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Abstract 



Pebble bed reactors (PBR) have moving graphite fuel pebbles. This unique 
feature provides advantages, but also means that simulation of the reactor 
requires understanding the typical motion and location of the granular flow 
of pebbles. 

This dissertation presents a method for simulation of motion of the peb- 
bles in a PBR. A new mechanical motion simulator, PEBBLES, efficiently 
simulates the key elements of motion of the pebbles in a PBR. This model 
simulates gravitational force and contact forces including kinetic and true 
static friction. It's used for a variety of tasks including simulation of the 
effect of earthquakes on a PBR, calculation of packing fractions, Dancoff fac- 
tors, pebble wear and the pebble force on the walls. The simulator includes 
a new differential static friction model for the varied geometries of PBRs. A 
new static friction benchmark was devised via analytically solving the me- 
chanics equations to determine the minimum pebble-to-pebble friction and 
pebble-to-surface friction for a five pebble pyramid. This pyramid check as 
well as a comparison to the Janssen formula was used to test the new static 
friction equations. 

Because larger pebble bed simulations involve hundreds of thousands of 
pebbles and long periods of time, the PEBBLES code has been parallehzed. 
PEBBLES runs on shared memory architectures and distributed memory 
architectures. For the shared memory architecture, the code uses a new 0(n) 
lock-less parallel collision detection algorithm to determine which pebbles are 
likely to be in contact. The new collision detection algorithm improves on 
the traditional non-parallel 0(n log(n)) collision detection algorithm. These 
features combine to form a fast parallel pebble motion simulation. 

The PEBBLES code provides new capabilities for understanding and op- 
timizing PBRs. The PEBBLES code has provided the pebble motion data 
required to calculate the motion of pebbles during a simulated earthquake. 



xii 



xiii 



The PEBBLES code provides the abihty to determine the contact forces and 
the lengths of motion in contact. This information combined with the proper 
wear coefficients can be used to determine the dust production from mechan- 
ical wear. These new capabihties enhance the understanding of PBRs, and 
the capabilities of the code will allow future improvements in understanding. 



ABSTRACT 



Chapter 1 
Introduction 



1.1 Pebble Bed Reactors Introduction 

Pebble bed nuclear reactors are a unique reactor type that have been pro- 
posed and used experimentally. Pebble bed reactors were initially developed 
in Germany in the 1960s when the AVR demonstration reactor was created. 
The 10 megawatt HTR-10 reactor achieved first criticality in 2000 in China 
and future reactors are planned. In South Africa, Pebble Bed Modular Reac- 
tor Pty. Ltd. was designing a full scale pebble bed reactor to produce process 
heat or electricity. 

Pebble bed nuclear reactors use graphite spheres (usually about 6 cm in 
diameter) for containing the fuel of the reactor. The graphite spheres encase 
smaller spheres of TRistructural-ISOtropic (TRISO) particle fuel. Unlike 
most reactors, the fuel is not placed in an orderly static arrangement. Instead, 
the graphite spheres are dropped into the top of the reactor, travel randomly 
down through the reactor core, and are removed from the bottom of the 
reactor. The pebbles are then possibly recirculated depending on the amount 
of burnup of the pebble and the reactor's method of operation. 

The first pebble bed reactor was conceived in 1950s in the West Germany 
using helium gas-cooling and spherical graphite fuel elements. Construc- 
tion on the Arbeitsgemeinschaft Versuchsreaktor (AVR) 15 MWe reactor 
was started in 1959 at the KFA Research Centre Jiilich. It started opera- 
tion in 1967 and continued operation for 21 years until 1988. The reactor 
operated with an outlet temperature of 950°C. The AVR demonstrated the 
potential for the pebble bed reactor concept. Over the course of its operation. 
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loss-of-coolant experiments were successfully performed. 

The second pebble bed reactor was the Thorium High Temperature Reac- 
tor (THTR). This reactor was built in West Germany for an electric utility. 
It was a 300 MWe plant that achieved full power in September 1986. In Oc- 
tober 1988, when the reactor was shutdown for maintenance, 35 bolt heads 
were found in the hot gas ducts leading to the steam generators. The deter- 
mination was made that the plant could be restarted, but funding difficulties 



prevented this from occurring and the reactor was decommissioned (Good- 



john 1991) 



The third pebble bed reactor to be constructed and the only one that is 
currently operable is the 10 MWt High Temperature Reactor (HTR-10). This 
reactor is at the Tsinghua University in China. Construction was started in 
1994 and reached first criticality in December 2000. This reactor is helium 



cooled and has an outlet temperature of 700°C (Wu et al. , 2002; Xu and Zuo 



2002) 



The use of high temperature helium cooled graphite moderated reactors 
with TRISO fuel particles have a number of advantages. A TRISO parti- 
cle consists of spherical fuel kernel (such as uranium oxide) surrounded by 
four concentric layers: 1) a porous carbon buffer layer to accommodate fis- 
sion product gases which limits pressure on the outer layers, 2) an interior 
pyrolytic carbon layer, 3) a layer of silicon carbide, and 4) an outer layer of 
pyrolytic carbon. The pyrolytic layers shrink and creep with irradiation, par- 
tially offsetting the pressure from the fission products in the interior as well 
as helping contain the fission gases. The silicon carbide acts as a containment 



mechanism for the metallic fission products. (Miller et al. , 2002) These layers 
provide an in-core containment structure for the radioactive fuel and fission 
products. 

The high temperature gas reactors have some advantages over conven- 
tional light water reactors. First, the higher outlet temperatures allow higher 
Carnot efficiency to be obtainecQ Second, the higher temperatures can be 
used for process heat, which can reduce the use of methane. Third, the 
high temperature under which TRISO particles can operate allows for the 
exploitation of the negative temperature coefficient to safely shutdown the re- 
actor without use of control rods|^ Fourth, the higher temperature is above 

^Thc outlet temperatures for pebble bed reactors have ranged from 700 °C to 950 °C, 
compared to typical outlet temperatures on the order of 300° C for light water reactors, so 
the intrinsic Carnot efficiency is higher. 

^Control rods are needed for a cold shutdown, however. 
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the annealing temperature for graphite, which safely removes Wigner en- 
ergjj^ These are advantages of both prismatic and pebble bed high temper- 



ature reactors. (Gougar et al. , 2004; Wu et al. , 2002) 



Pebble bed reactors, unlike most other reactors types, have moving fuel. 
This provides advantages but complicates modeling the reactors. A key ad- 
vantage is that pebble bed reactors can be refueled online, that is, reactor 
shutdown is not needed for refueling. As a consequence, the reactors have 
low excess reactivity, as new pebbles can be added or excess pebbles removed 
to maintain the reactor at critical. The low excess reactivity removes the 
need for burnable poisons. A final advantage is that the moving fuel allows 
the pebble bed to be run with optimal moderation, where both increases 
and decreases in the fuel-to-moderator ratio cause reduction in reactivity. 



Ougouag et al. (2004) discuss the advantages of optimal moderation includ- 
ing improved fuel utilization. However, because the fuel is moving, many 
traditional methods of modeling nuclear reactors are inapplicable without a 
method for quantifying the motion. Hence, there is a need for development 
of methods usable for pebble bed reactor modeling. 



1.2 Dissertation Introduction 



This dissertation describes a computer code, PEBBLES, that is designed to 
provide a method of simulating the motion of the pebbles in a pebble bed 
reactor. 

Chapter |4] provides the details of how the simulation works. Chapter [5] 
has a new static friction model developed for this dissertation. 

Several checks have been made of the code. Figure 3.1 compares the 



PEBBLES simulation to experimentally determined radial packing fractions. 



Section 5A_ describes a new analytical benchmark that was used to test the 
static friction model in PEBBLES. Section 15.21 uses the Janssen model to 
test the static friction in a cylindrical vat. 

Motivating all the above are the new applications, including Dancoff fac- 



tors (8. LI), calculating the angle of repose (8.1.2) and modeling an earth- 



quake in section 8.2 



■^The accumulation of Wigner energy led to the Windscale fire in that lower temperature 
graphite reactor. 



Chapter 2 
Motivation 



Most nuclear reactors have fixed fuel, including typical light water reactors. 
Some reactor designs, such as non-fixed fuel molten salt reactors, have fuel 
that is in fluid flow. Most designs for pebble bed reactors, however, have 
moving granular fuel. Since this fuel is neither fixed nor easily treatable as a 
fluid, predicting the behavior of the reactor requires the ability to understand 
the characteristics of the positions and motion of the pebbles. For example, 
predicting the probability of a neutron leaving one TRISO's fueled region and 
entering another fueled region depends on the typical locations of the pebbles. 
A second example is predicting the effect of an earthquake on the reactivity 
of the pebble bed reactor. This requires knowing how the positions of the 
pebbles in the reactor change from the forces of the earthquake. Accurate 
prediction of the typical features of the flow and arrangement of the pebbles in 
the pebble bed reactor would be highly useful for their design and operation. 

The challenge is to gain the ability to predict the pebble flow and pebble 
positions for start-up, steady state and transient pebble bed reactor opera- 
tion. 

The objective of the research presented in this dissertation is to provide 
this predicting ability. The approach used is to create a distinct element 
method computer simulation. The simulation determines the locations and 
velocities of all the pebbles in a pebble bed reactor and can calculate needed 
tallies from this data. Over the course of creating this simulation, various 
applications of the simulation were performed. These models allow the op- 
eration of the pebble bed reactor to be better understood. 
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Chapter 3 
Previous work 



Because the purpose of this dissertation is to produce a high fidehty simula- 
tion that can provide predictions of the pattern and flow of pebbles, previous 
efforts to simulate granular methods and packing were examined. A variety 
of simulations of the motion of discrete elements have been created for differ- 
ent purposes. Lu et al. (2001) applied a discrete element method (DEM) to 
determine the characteristics of packed beds used as fusion reactor blankets. 



JuUien et al. ( 1992 ) used a DEM to determine packing fractions for spheres 



using different non-motion methods. 



Soppe (1990) used a rain method to 
The rain method ran- 



determine pore structures in different sized spheres 
domly chooses a horizontal position, and then lowers a sphere down until it 
reaches other existing spheres. This is then repeated to fill up the container. 



Freund et al. (2003) used a rain method for fluid flow in chemical processing. 



The use of non-motion pebble packing methods provide an approximation 
of the positions of the pebble. Unfortunately, non-motion methods will tend 
to either under pack or over pack (sometimes both in the same model). 
For large pebble bed reactors, the approximately ten-meter height of the 
reactor core will result in different forces at the bottom than at the top. 
This will change the packing fractions between the top and the bottom, so 
without key physics, including static friction and the transmittal of force, 
non-motion physics models will not even be able to get correct positional 
information. Non-physics based modeling can not be used for predicting the 
effect of changes in static friction or pebble loading methods even if only the 
position data is required. 

The initial PEBBLES code for calculation of pebble positions minimized 
the sum of the gravitational and Hookes' law potential energies by adjusting 
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pebble positions. However, that simulation was insufficient for determining 
flow and motion parameters and simulation of earthquake packing. 

Additional references addressing full particle motion simulation were eval- 
uated. Kohring ( 1995 ) created a 3-D discrete element method simulation to 
study diffusional mixing and provided detailed information on calculating 
the kinetic forces for the simulation. The author describes a simple method 
of calculating static friction. Haile (1997) discusses both how to simulate 
hard spheres and soft spheres using only potential energy. The soft sphere 
method in Haile proved useful for determining plausible pebble positions, but 
is insufficient for modeling the motion. Hard spheres are simulated by calcu- 
lating the collision results from conservation laws. Soft spheres are simulated 
by allowing small overlaps, and then having a resulting force dependent on 
the overlap. Soft spheres are similar to what physically happens, in that the 
contact area distorts, allowing distant points to approach closer than would 
be possible if the spheres were truly infinitely hard and only touched at one 
infinitesimal point. Hard spheres are impractical for a pebble bed due to 
the frequent and continuous contact between spheres so soft spheres are used 
instead. The dissertation by Ristow (1998) describes multiple methods for 
simulation of granular materials. On Ristow's list of methods was a model 
similar to that used as the kernel of the work supporting this dissertation. 
Ristow's dissertation mentioned static friction and provided useful references 
that will be discussed in Section [3^ 

To determine particle fiows. Wait (2001) developed a discrete element 
method that included only dynamic friction. Concurrently with this dis- 
sertation research, Rycroft et al. (2006b) used a discrete element method, 
created for other purposes, to simulate the fiow of pebbles through a pebble 
bed reactor. 

Multiple other discrete element codes have been created and PEBBLES 
is similar to several of the full motion models. For most of the applications 
discussed in this dissertation, only a model that simulates the physics with 
high fidelity is useful. The PEBBLES dynamic friction model is similar to 
the model used by Wait or Rycroft, but the static friction model incorporates 
some new improvements that will be discussed later. 

In addition to simulation by computer, other methods of determining the 
properties of granular fiuids have been used. Bedenig et al. (1968) used a 
scale model to experimentally determine residence spectra (the amount of 
time that pebbles from a given group take to pass through a reactor) for 
different exit cone angles. Kadak and Bazant (2004) used scale models and 
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small spheres to estimate the flow of pebbles through a full scale pebble 
bed reactor. These researchers also examined the mixing that would occur 
between different radial zones as the pebbles traveled downward. Bernal 



et al. (1960) carefully lowered steel spheres into cylinders and shook the 



cylinders to determine both loose and dense packing fractions. The packing 
fraction and boundary density fluctuations were experimentally measured by 
Benenati and Brosilow (1962). The Benenati and Brosilow data have been 



used to verify that the PEBBLES code was producing correct boundary 
density fluctuations (See Figure 3.1). Many experiments were performed 
in the designing and operating of the AVR reactor to determine relevant 



properties such as residence times and optimal chute parameters (Baumer 



et al. , 1990). These experiments provide data for testing the implementation 



of any computational model of pebble flow. 

The PEBBLES simulation uses elements from a number of sources and 
uses standard classical mechanics for calculating the motion of the pebbles 
based on the forces calculated. The features in PEBBLES have been chosen 
to implement the necessary fldelity required while allowing run times small 
enough to accommodate hundreds of thousands of pebbles. The next sections 
will discuss handling static friction. 



3.1 Static Friction Overview 

Static friction is an important effect in the movement of pebbles and their 
locations in pebble bed reactors. This section briefly reviews static friction 
and its effects in pebble bed reactors. Static friction is a force between two 
contacting bodies that counteracts relative motion between them when they 
are moving sufficiently slowly (Marion and Thornton, 2004). Macroscopically, 



the maximum magnitude of the force is proportional to the normal force with 
the following equation: 



\Fs\ </i|Fj 



(3.1) 

is the static friction force 



where /i is the coefficient of static friction, 
and F± is the normal (load) force. 

Static friction results in several effects on granular materials. Without 
static friction, the angle of the slope of a pile of a material (angle of repose) 



would be zero(Duran, 1999). Static friction also allows 'bridges' or arches to 



be formed near the outlet chute. If the outlet chute is too small, the bridging 
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Figure 3.1: Comparison Between PEBBLES Outputs and Benenati 
and Brosilow Data 



will be stable enough to clog the chute. Static friction will also transfer force 
from the pebbles to the walls. This will result in lower pressure on the walls 



than would occur without static friction(Sperl, 2006 Walker, 1966). 



For an elastic sphere, static friction's counteracting force is the result of 
elastic displacement of the contact point. Without static friction, the contact 
point would slide as a result of relative motion at the surface. With static 
friction, the spheres will experience local shear that distorts their shape so 
that the contact point remains constant. This change will be called stuck- 
slip, and continues until the counteracting force exceeds fiF±. When the 
counteracting force exceeds that value, the contact point changes and slide 
occurs. The mechanics of this force with elastic spheres were investigated by 
Mindlin and Deresiewicz (1953). Their work created exact formulas for the 
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force as a function of the past relative motion and force. 



3.2 Simulation of Mechanics of Granular Ma- 
terial 



Many simulations of granular materials incorporating static friction have 
been devised. Cundall and Strack (1979) developed an early distinct ele- 



ment simulation of granular materials that incorporated a computationally 
efficient static friction approximation. Their method involved integration of 
the relative velocity at the contact point and using the sum as a proxy for 
the current static friction force. Since their method was used for simulation 
of 2-D circles, adaptation was required for 3-D granular materials. One key 
aspect of adaptation is determining how the stuck-slip direction changes as 
a result of contacting objects' changing orientation. 



Vu-Quoc and Zhang (1999) and Vu-Quoc et al. (2000) developed a 3-D 



discrete-element method for granular flows. This model was used for simu- 
lation of particle flow in chutes. They used a simplification of the Mindlin 
and Deresiewicz model for calculating the stuck-slip magnitude, and project 
the stuck-slip onto the tangent plane each time-step to rotate the stuck-slip 
force direction. This correctly rotates the stuck-slip, but requires that the 
rotation of the stuck-slip be done as a separate step since it not written in a 
differential form. 



Silbert et al. (2001) and Landry et al. (2003) describe a 3-D differential 



version of the Cundall and Strack method. The literature states that particle 
wall interactions are done identically. The amount of computation of the 
model is less than the Vu-Quoc, Zhang and Walton model. This model 
was used for modeling pebble bed flow(Rycroft et al. 2006a|[b ). This model 



however, does not specify how to apply their differential version to modeling 
curved walls. 



Chapter 4 
Mechanics Model 



The PEBBLES simulation calculates the forces on each individual pebble. 
These forces are then used to calculate the subsequent motion and position 
of the pebbles. 



4.1 Overview of Model 

The PEBBLES simulation tracks each individual pebble's velocity, position, 
angular velocity and static friction loadings. The following classical mechan- 
ics differential equations are used for calculating the time derivatives of those 
variables: 



dt rrii 
dpi 

^^^^ 
duji _ J2i^j F\iij X nuij + F|| 



dt L 



(4.1) 
(4.2) 
(4.3) 



= S(Fxij, Vi, Vj, Pi, pj, Sij) (4.4) 

where Fj^ is the force from pebble j on pebble i, Yd is the force of the 
container on pebble g is the gravitational acceleration constant, is the 
mass of pebble i, is the velocity of pebble i, pi is the position vector for 
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pebble i, ui is the angular velocity of pebble i, F||jj is the tangential force 
between pebbles i and j, F±ij is the perpendicular force between pebbles i 
and j, Ti is the radius of pebble i, li is the moment of inertia for pebble i, 
F||ci is the tangential force of the container on pebble i, Ad is the unit vector 
normal to the container wall on pebble i, fiij is the unit vector pointing from 
the position of pebble i to that of pebble j, Sij is the current static friction 
loading between pebbles i and j, and S is the function to compute the change 
in the static friction loading. The static friction model contributes to the F||jj 
term which is also part of the Fjj term. Figure 4.1 illustrates the principal 
vectors with pebble i going in the Vj direction and rotating around the Ui 
axis, and pebble j going in the Vj direction and rotating around the uj axis. 




Figure 4.1: Principle Vectors in the Interaction of Two Pebbles 

The mass and moment of inertia are calculated assuming spherical sym- 
metry with the equations: 

4 

m = -7i [p,rl + po{rl - r^)] (4.5) 

/ = ^TT [p^rl + po{rl - rl)] (4.6) 

where Tc is the radius of inner (fueled) zone of the pebble, is the radius 
of whole pebble, pc is the average density of center fueled region and po is 
the average density of outer non-fueled region. 
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The dynamic (or kinetic) friction model is based on the model described 



by Wait (2001). Wait's and PEBBLES model calculate the dynamic fric- 



tion using a combination of the relative velocities and pressure between the 



pebbles, as shown in Equations (4.7) and (4.8): 

F±ij = hlijhij - C±\±ij, lij > (4.7) 

Fd\\ij = -mm(/i|F_Lij|,C|||v||ij|)v||ij-,/ij > (4.8) 

where Cy is the tangential dash-pot constant, C± is the normal dash- 
pot constant, F±ij is the normal force between pebbles i and j, Fd\\ij is the 
tangential dynamic friction force between pebbles i and j, h is the normal 
Hooke's law constant, l^j is the overlap between pebbles i and j, vyjj is 
the component of the velocity between two pebbles perpendicular to the 
line joining their centers, is the component of the velocity between two 
pebbles parallel to the line joining their centers, Vjj is the relative velocity 
between pebbles i and j and /i is the kinetic friction coefficient. Equations 



(4.9 4.12) relate supplemental variables to the primary variables: 

Fij = F^,j + F||,j (4.9) 

v±ij = i^ij ■ n^j)nij (4.10) 

= v^j' - v^ij (4.11) 

Vjj = (vj + UJiX TiUij) - (vj + iOj X rjhji) (4.12) 

The friction force is then bounded by the friction coefficient and the 
normal force, to prevent it from being too great: 

F/lHi = Fs\\ij + Fd\\ij (4.13) 

F|,i, = mm(/i|F^,j|, \F f\\ij\)F f\\,j (4.14) 

where Fs||ij is the static friction force between pebbles i and j, Fd\\ij is 
the kinetic friction force between pebbles i and j, hg is the coefficient for 
force from slip, Sij is the slip distance perpendicular to the normal force 



4.2. INTEGRATION 
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between pebbles i and j, Vmax is the maximum velocity under which static 
friction is allowed to operate, and fi is the static friction coefficient when the 
velocity is less than v^a^ and the kinetic friction coefficient when the velocity 
is greater. These equations fully enforces the first requirement of a static 
friction method, |Fs| < /i|Fj^|. 

4.2 Integration 

When all the position, linear velocity, angular velocity and slips are com- 
bined into a vector y, the whole computation can be written as a differential 
formulation in the form: 

y' = f(t,y) (4.15) 
y{to) = yo (4.16) 

This can be solved by a variety of methods with the simplest being Euler's 
method: 

yi = yo + Aif(i,yo) (4.17) 

In addition, both the Runge-Kutta method and the Adams-Moulton 
method can be used for solving this equation. These methods improve the 
accuracy of the simulation. However, they do not improve the wall-clock 
time at the lowest stable simulation, since the additional time required for 
computation negates the advantage of being able to use somewhat longer 
time-steps. In addition, when running on a cluster, more data needs to be 
transferred since the methods allow non-contacting pebbles to affect each 
other in a single 'time-step calculation'. 

4.3 Geometry Modeling 

For any geometry interaction, two things need to be calculated, the overlap 
distance I (or, technically, the mutual approach of distant points) and the 
normal to the surface n. The input is the radius of the pebble r and the 
position of the pebble, p with components p^, Py, and Pz 
For the floor contact this is: 
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I — {pz — r) — floor Jocation 
n = z 

For cylinder contact on the inside of a cylinder this is: 
I — {pr + r) — cylinder ^radius 

-Px^ , -Py ^ 

n — X H -y 

pr pr 

For cyhnder contact on the outside of a cyhnder this is: 

^JpI + pI 

I = cylinder .radius + r — pr 

Px. , Py. 

n = — X H — -y 
pr pr 

For contact on the inside of a cone defined by the radius — mz + h: 



pr=^Jpl + Pl (4.26) 

z^ = ^^P^-^)^' (4.27) 
+ 1 

Tc = mzc + h (4.28) 

Xc = {rc/pr)Py (4.29) 

Vc = {rc/pr)p^ (4.30) 

XcX + ydj + ZcZ (4.31) 

d = p - c (4.32) 

Z = |d| +r,rc <pr (4.33) 

n — —d, Tc < pr (4.34) 

I — r — \d\,rc >= pr (4.35) 

n = d,rc >= pr (4.36) 



(4.18) 
(4.19) 



(4.20) 
(4.21) 
(4.22) 



(4.23) 
(4.24) 
(4.25) 
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These equations are derived from minimizing the distance between the 
contact point c and the pebble position p. 

For contact on a plane defined hj ax + by + cz + d = where the equation 
has been normalized so that + 6^ + = 1, the following is used: 



dp = apx + bpy + cpz + d (4.37) 
l = r-dp (4.38) 
h = ax + by + cz (4.39) 

Combinatorial geometry operations can be done. Intersections and unions 
of multiple geometry types are done by calculating the overlaps and normals 
for all the geometry objects in the intersection or union. For an intersection, 
where there is overlap on all the geometry objects, then the smallest overlap 
and associated normal are kept, which may be no overlap. For a union, the 
largest overlap and its associated normal are kept. 

For testing that a geometry is correct, a simple check is to fill up the 



geometry with pebbles using one of the methods described in Section AA_ 
and then make sure that linear and angular energy dissipate. Many geome- 
try errors will show up by artificially creating extra linear momentum. For 
example, if a plane is only defined at the top, but it is possible for pebbles to 
leak deep into the bottom of the plane, they will go from having no overlap 
to a very high overlap, which will give the pebble a large force. This results 
in extra energy being added each time a pebble encounters the poorly defined 
plane, which will show up in energy tallies. 



4.4 Packing Methods 

The pebbles are packed using three main methods. The simplest creates a 
very loose packing with an approximately 0.15 packing fraction by randomly 
choosing locations, and removing the overlapping ones. These pebbles then 
allowed to fall down to compact to a realistic packing fraction. 



The second is the PRIMe method developed by Kloosterman and Ougouag 



(2005). In this method large numbers of random positions (on the order of 



100,000 more than will fit) are generated. The random positions are sorted 



by height, and starting at the bottom, the ones that fit are kept. Figure 4.2 



illustrates this process. This generates packing fractions of approximately 
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0.55. Then they are allowed to fall to compact. This compaction takes less 
time than starting with a 0.15 packing fraction. 



A"''i-' -4 yr: V—..—-. 




Figure 4.2: PRIMe Method Illustration 



The last method is to automatically generates virtual chutes above the 
bed where the actual inlet chutes are, and then loads the pebbles into the 



chutes. Figure 4^ shows this in progress. This allows locations that have 
piles where the inlet chutes are, but can be done quicker than a recirculation. 
The other two methods generate fiat surfaces at the top, which is unrealistic, 
since the surface of a recirculated bed will have cones under each inlet chute. 



4.5 Typical Parameters 

The typical parameters used with the PEBBLES code are described in Table 
14.11 Alternative numbers will be described when used. 
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Virtual 




Figure 4.3: Virtual Chute Method 



Table 4.1: Typical Constants used in Simulation 



Constant 


Value 


Gravitational Acceleration g 


9.8 m/s^ 


radius of pebbles r 


0.03 m 


Mass of Pebble m 


0.2071 kg 


Moment of Inertia / 


7.367e-05 kg m^ 


Hookc's law constant h 


1.0e6 


Dash-pot constants Cy and C± 


200.0 


Kinetic Friction Coefficient n 


0.4 or sometimes 0.25 


Static Friction Coefficient fj,s 


0.65 or sometimes 0.35 


Maximum static friction velocity v^ax 


0.1 m/s 



Chapter 5 

A New Static Friction Model 



The static friction model in PEBBLES is used to calculate the force and mag- 
nitude of the static friction force. Other models have been created before to 
calculate static friction, but the PEBBLES model provides the combination 
of being a differential model (as opposed to one where the force is rotated as 
a separate step) and being able to handle the type of geometries that exist 
in pebble bed reactors. 

The static friction model has two key requirements. First, the force from 
stuck-slip must be updated based on relative motion of the pebbles. Second, 
the current direction of the force must be calculated since the pebbles can 
rotate in space. 



5.0.1 Use of Parallel Velocity for Slip Updating 



For elastic spheres, the true method of updating the stuck-shp force is to use 



the method of Mindlin and Deresiewicz (1953). This method requires com- 



putationally and memory intensive calculations to track the forces. Instead, 
a simpler method is used to approximate the force. This method, described 



by Cundall and Strack (1979) uses the integration of the parallel relative ve- 



locity as the displacement. The essential idea is that the farther the pebbles 
have stuck-slipped at the contact point, the greater the counteracting static 
friction force needs to be. This is what happens under more accurate models 
such as Mindlin and Deresiewicz. There are two approximations imposed 
by this assumption. First, the amount the force changes is independent of 
the normal force. Second, the true hysteretic effects that are dependent on 
details of the loading history are ignored. For simulations where the exact 
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dynamics of static friction are important, these could potentially be serious 
errors. However, since static friction only occurs when the relative speed 
is low, the dynamics of the simulation usually are unimportant. Thus, for 
most circumstances, the following approximation can be used for the rate of 
change of the magnitude and non-rotational change of the stuck-slip: 



5.0.2 Adjustment of Oversize Slips 

The slips can build up to unrealistically large amounts, that is, greater than 



/i|Fx|; equation 5.1 places no limit on the maximum size of the slip. The 
excessive slip is solved at two different locations. First, when the frictions 
are added together to determine the total friction they are limited by /i|Fx| 



in equation (4.14). This by itself is insufficient, because the slip is storing 
potential energy that appears anytime the normal force increases. This man- 
ifests itself by causing vibration of the pebbles to continue for long periods of 
time. Two methods for fixing the hidden slip problem are available in PEB- 
BLES. The simplest drops any slip that exceeds the static friction threshold 
(or an input parameter value somewhat above the static friction threshold 
so small vibrations do not cause the slip to disappear). 

The second method used in PEBBLES is to decrease the slip that is over 
a threshold value. If the slip is too great, a derivative that is the opposite 
as the current slip is added as an additional term in the slip time derivative. 
This occurs in the following additional term: 

= -R{\sij\ - SsdfJ'\F±ij\)s^ij (5.2) 

In this R(a;) is the ramp function (which is a; if x > and otherwise), 
Ssd is a constant to select how much the slip is allowed to exceed the static 
friction threshold (usually 1.1 in PEBBLES). This derivative adder is used 
in most PEBBLES runs since it does allow vibrational energy to decrease, 
yet does not cause the pyramid benchmark to fail like complete removal of 
too great slips does. 

When using non-Euler integration methods, the change in slip is calcu- 
lated multiple times. Each time it is calculated, it might be set to be zeroed. 
In the PEBBLES code, if any of the added up slips for a given contact were 
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set to be zeroed, the final slip is zeroed. This is not an ideal method, but it 
works well enough. 



5.0.3 Rotation of Stuck-Slip 

The static friction force must also be rotated so that it is in the plane of 
contact between the two pebbles. When there is a difference between the 
pebbles' center velocities, which changes in the relative pebble center loca- 
tion, change in the direction in the stuck-slip occurs. That is: 

Pm+l — Pjn+l ~ Pin — Pjn + (Vjn — Vj,-„)At (5-3) 

First, let njj„ = Pi — pj and duijn = Vj — Vj. The cross product —duijn x 
Uijn is perpendicular to both n and dn and signed to create the axis around 
which s is rotated in a right-handed direction. Then, using the cross product 
of the axis and s, — (rfrijj x njj„) x s^n gives the correct direction that s 
should be increased. 

Next, determine the factors required to make the differential the proper 
length. By cross product laws. 




Figure 5.1: Static Friction Vectors 



The goal is to rotate s by angle af which is the 'projection' into the 
proper plane of the angle a that n rotates by. Since the direction has been 
determined, for simplicity the figure leaves the indexes off, and concentrates 
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on determining the lengths. In Figure 5A^. s is the old slip vector, s/ is the 
new slip vector, n is the vector pointing from one pebble to another. The 
vector dnAt is added to n to get the new n/, n + dnAt. The initial condition 
is that s and n are perpendicular. The final conditions are that s/ and n/ 
are perpendicular, and that s and s/ are the same length and that s/ is the 
closest vector to s as it can be while satisfying the other conditions. There 
is no requirement that s or s/ are coplanar with dnAt (otherwise a/ would 
be equal to a). From the law of sines we have: 

\dnAt\ |n| 

(5.5) 



sin a sin 6 
so 



|cinAt|sin^ 

sma = 1 — I (5.0) 

n 



In Figure 5.2 the projection to the correct plane occurs. First by using 
the length of s is projected to the plane. Note that is the angle both to 
s and to s/. So, the length of the line on the dn x n plane is |s| sin0, and 
the length of the straight line at the end of the triangle is |s| sine/) sin a (note 
that the chord length is |s|(sin0)Q;, but as At approaches the other can 
be used). From these calculations, the length of the dsAt can be calculated 
with the following formula: 

, . Isl sin 0|(inAt| sin 

dsAt = ^ — ^ 5.7 

|n| 

Since | — {diiij x njj„) x Sjj„| = |(injj| |njj„| |Sjj„| sin 6' sine/) the formula for 
the rotation is: 

As a differential equation this is: 



rfSij _ [((Vi - Vj) X (pi - pj-)) X s. 



dt |Pi-PjP 
By the vector property a x (6 x c) = h{a-c)—c{a-\)) and since (pj — pj)-Sjj 



(5.9) 



0, this can be rewritten as the version in Silbert et al. (2001): 



^ ^ (P»-Pj)(Su-(vi-Vj)) 
dt |Pi-PjP 



(5.10) 
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|s|sin CD sin a 



Figure 5.2: Projections to ds 



5.0.4 Differential Equation for Surface Slip Rotating 

It might seem that the wall interaction could be modeled the same way as the 
pebble-to-pebble interaction. For sufficiently simple wall geometries this may 
be possible, but actual pebble bed reactor geometries are more complicated, 
and violate some of the assumptions that underpin the derivation. For a flat 
surface, there is no rotation, so the formula can be entirely dropped. For a 
spherical surface, it would be possible to measure the curvature at pebble to 
surface contact point in the direction of relative velocity to the surface. This 
curvature could then be used as an effective radius in the pebble-to-pebble 
formulas. 

The pebble reactor walls have additional features that violate assumptions 
made for the derivation. For surfaces such as a cone, the curvature is not in 
general constant, because the path can follow elliptical curves. As well, the 
curvature has discontinuities where different parts of the reactor join together 
(for example, the transition from the outlet cone to the outlet chute). At 
these transitions, the assumption that the slip stays parallel to the surface 
fails because the slip is parallel to the old surface, but the new surface has a 
different normal. 

Because of the complications with using the pebble to pebble interaction, 
PEBBLES uses an approximation of the "rotation delta." This is similar 



to the Vu-Quoc and Zhang (1999) method of adjusting the slip so that it is 



parallel to the surface each time. Each time when the slip is used, a tempo- 
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rary version of the slip that is properly aligned to the surface is computed 
and used for calculating the force. As well, a rotation to move the slip more 
parallel to the surface is also computed. 

The rotation is computed as follows. Let the normal direction of the wall 
at the point of contact of the pebble be n, and the old stuck-slip be s. Let a 
be the angle between n and s. n x s is perpendicular to both n and s and so 
this cross product is the axis that needs to be rotated around. Then (n x s) x s 
is perpendicular to this vector, so it is either pointing directly towards n if 
a is acute or directly away from n if a is obtuse. To obtain the correct 
direction, this vector is multiplied by the scalar s ■ n which has the correct 
sign from cos a. The magnitude of (s ■ n)[(n x s) x s] needs to be determined 
for reasonableness. We define the angle 6, which is between (n x s) and s. 
By these definitions the magnitude is (|s||n| cosa)[(|n||s| sina)|s| sin6]. h is 
a right angle since n x s is perpendicular to s, so sin 6 = 1. Collecting terms 
gives the magnitude as |sp|np cosasina which is divided by |n x s||n||s| to 
give the full term the magnitude |s| cos a. This is the length of the vector that 
goes from s to the plane perpendicular to n. This produces equation |5.11 



which can be used to ensure that the wall stuck-slip vector rotates towards 
the correct direction. 



ds 
dt 



s ■ n 



[(n X s) X s] 
In X s| Inl Isl 



(5.11) 




old s 



Figure 5.3: Static Friction Vectors for Wall 
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5.1 Testing of Static Friction Model With Pyra- 
mid Test 

Static friction is an important physical feature in the implementation of me- 
chanical models of pebbles motion in a pebble bed, and checking its correct- 
ness is important. A pyramid static friction test model was devised as a 
simple tool for verifying the implementation of a static friction model within 
the code. The main advantages of the pyramid test are that the model test 
is realistic and that it can be modeled analytically, providing an exact basis 
for the comparison. The test benchmark consists of a pyramid of five spheres 
on a flat surface. This configuration is used because the forces acting on 
each pebble can be calculated simply and the physical behavior of a model 
with only kinetic friction is fully predictable on physical and mathemati- 
cal grounds: with only kinetic friction and no static friction, the pyramid 
will quickly flatten. Even insufficient static friction will result in the same 
outcome. The four bottom spheres are arranged as closely as possible in a 



square, and the fifth sphere is placed on top of them as shown in Fig. |5.4 

Top View 

Side View 





Figure 5.4: Sphere Location Diagram 

The lines connecting the centers of the spheres form a pyramid with 



sides 2R, as shown in Fig. 5.5, where R is the radius of the spheres. The 



length of a in the figure is and because b is part of a right triangle, 
{2Rf - i^f = = 4i?2 _ 4|! = 2R^, so 6 has the same length as a, 
and thus the elevation angle for all vertexes of the pyramid are 45° from 
horizontal. 

Taking for origin of the coordinates system the projection of the pyramid 
summit onto the ground, the locations (coordinates) of the sphere centers 



are given in Table 5.1 
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Figure 5.5: Pyramid Diagram 
Table 5.1: Sphere location table. 



X 


Y 


Z 


-R 


-R 


R 


R 


-R 


R 


-R 


R 


R 


R 


R 


R 








R{1 + V2) 



5.1.1 Derivation of Minimum Static Frictions 

The initial forces on the base sphere are the force of gravity mg, and the 



normal forces Tn and Fn as shown in Fig. 5^ This causes initial stuck-slip 
which will cause Fs to develop to counter the slip, and Ts to counter the 
rotation of the base sphere relative to the top sphere. The top sphere will 
have no rotation because the forces from the four spheres will be symmetric 
and counteract each other. 

The forces on the base sphere are: 

Tn - Normal force from the top sphere 

Ts - Static friction force from the top sphere 

mg - Force of gravity on the base sphere 

Fn - Normal force from floor 

Fs - Static friction force from the floor 



Note that Fn is larger than Tn since Tn is only a portion of the mg 
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Figure 5.6: Force Diagram 



force since the top sphere transmits (and sphts) its force onto all four base 
spheres. 

There are three requirements for a base sphere to be non- accelerated. 
If a base sphere is not rotating than there is no torque, so: 

|Fs| = |Ts| (5.12) 

The resultant of all forces must also be zero in the x and the y direc- 
tion (vector notation dropped since they are in one dimension and therefore 
scalars) as follows: 

-Fs - Tsx + Tnx = (5.13) 



—mg — Tsy — Tny + Fn = 



(5.14) 
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Since the angle of contact between a base sphere and the top sphere is 
45°, the following two equations hold (where Ts is the magnitude of Ts and 
Tn is the magnitude of Tn): 



Tsx = Tsy - 
Tnx = Tny 



Ts 
Tn 

71 



This changes equations 5.13 and 5.14 into: 



Ts Tn 
-Fs - ^ + ^ = C 

Ts Tn 

—mq — — + Fn 

^2 72 



Combining equation 5.12 and |5 . 1 7] provides : 



^ Ts Tn 
-Ts-^ + ^ 



Which gives the relation: 



Tn = Ts(V2 + 1) 



By the static friction Equation 3.1 



Ts ^ fXgpfiQ^^Tn 



(5.15) 
(5.16) 

(5.17) 
(5.18) 

(5.19) 

(5.20) 
(5.21) 



Combining equations 5.20 and 5.21 and simplifying gives the requirement 



that 

V2-I < fisphere (5.22) 

For use with testing, the static friction program can be tested twice with a 
sphere-to-sphere friction coefficient slightly above 0.41421... and one slightly 
below 0.41421.... In the first case the pyramid should be stable, and in the 
second case the top ball should fall to the fioor. 

Since j of the weight of the top pebble is on one of the base pebbles, the 
following holds: 



Fn = -mq 
4 



(5.23) 
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Combining 5.18 and 5.23 provides the following equation: 

mg Ts Tn 



4 ^2 ^2 







(5.24) 



Equations 5.17 and 5.24 can be added to produce 



-Fs - V2Ts + 



mg 



(5.25) 



Using 5.12 and 5.24 and solving for Fs gives the following value for Fs: 



Fs 



mg 



4(1 + V2) 



(5.26) 



By the static friction Equation 3.1 



Fs ^ ^surfaceF^- 

Substituting the values for Fs and Fn gives: 



mg 



4(1 + y/2) 



^ f^surface . mg 



5 

= 4' 



(5.27) 



(5.28) 



Simplifying provides the following relation for the surface-to-sphere static 
friction requirement: 



1 



5(1 + V2) 



< 



f^surfa 



(5.29) 



This can be used similarly to the other static friction requirement by 
setting the value slightly above 0.08284... and slightly below 0.08284... and 
making sure that it is stable with the higher value and not stable with the 
lower value. 

This test was inspired by an observation of lead cannon balls stacked into 
a pyramid. I tried to stack used glass marbles into a five ball pyramid and 
it was not stable. Since lead has a static friction coefficient around 0.9 and 
used glass has a much lower static friction, the physics of pyramid stability 
was further investigated, resulting in this benchmark test of static friction 
modeling. 



5.2. JANSSEN'S FORMULA COMPARISON 
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5.1.2 Use of Benchmark 

The benchmark test of two critical static friction coefficients is defined by the 
following equations. If both static friction coefficients are above the critical 
values, the spheres will form a stable pyramid. If either or both values are 
below the critical values the pyramid will collapse. 

IJ'criticalsurface = — 7^ ~ 0.08284 (5.30) 

5(1 + V2) 

^J'criticalsphere = — 1 ~ 0.41421 (5.31) 



To set up the test cases, the sphere locations from Table |5.1| should be 
used as the initial locations of the sphere. Then, static friction coefficients 
for the sphere-to-sphere contact and the sphere-to-surface contact are chosen. 
The code is then run until either the center sphere falls to the surface, or the 
pyramid obtains a stable state. There are three test cases that are run to 
test the model. 

!• f^surface f^criticalsurface ~t" ^ and ^sphere ^^ critical sphere ~t" ^ which should 

result in a stable pyramid. 

2. surf ace critical sur f ace ^ and ^sphere ^^criticalsphere ~t~ ^ which should 

result in a fall. 

3. surf ace f^criticalsurface ~t~ ^ and ^sphere f^criticalsphere ^ which should 

result in a fall. 

For soft sphere models, there are fundamental limits to how close the 
model's results can be to the critical coefficient. Essentially, as the critical 
coefficients are approached, the assumptions become less valid. For example, 
with soft (elastic) spheres, the force from the center sphere will distort the 
contact angle, so the actual critical value will be different. For the PEBBLES 
code, an e value of 0.001 is used. 



5.2 Testing of the Static Friction Model by 
Comparison with Janssen's Formula 



The pyramid static friction test is used as a simple test of the static friction 
model. Another test compares the static friction model against the Janssen 
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formula's behavior (Sperl 2006). This formula specifies the expected wall 



pressure as a function of depth. This formula only applies when the static 
friction is fully loaded, that is when F^l = /i|F_|_|. This condition is generally 



not satisfied until some recirculation has occurred. Figure 5.7 shows the 
normal force and the static friction force from a pebble to the wall. With the 
PEBBLES code, this is only satisfied after recirculation with lower values of 
the static friction coefficient /i. 




static 



Figure 5.7: Relevant Forces on Wall from Pebble 

The equation used to calculate the pressure on the region R from the 
normal force in PEBBLES is: 



P 



Rh2nr 



(5.32) 



in R 



where p is the pressure, Rh is the height of the region, and r is the radius 
of the cylinder. 

The equation for calculating the Janssen formula pressure is 



K — 2/ipp — 2/ippW fipp 



1 + 1 



fpg 

b2r 



P 



4/i 



wall 



1 — e 



(5.33) 
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where Upp is the pebble to pebble static friction coefficient, fiwaii is the 
pebble to wall, / is the packing fraction, p is the density, g is the gravitational 
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acceleration, and z is the depth that the pressure is being calculated. For 
the Figures 5^ and 5^, a packing fraction of 0.61 is used and a density 
of 1760 kg/m^ are used. There are 20,000 pebbles packed into a 0.5 meter 
radius cylinder, and 1,000 are recirculated before the pressure measurement 
is done. 

Figure 5^ compares the Janssen model with the PEBBLES simulation 
for static friction values of 0.05 and 0.15. For this case, the Janssen formula 
and the simulated pressures match closely. Figure 5^ compares these again. 
In this case, the 0.25 fi values only approximately match, and the 0.9 static 
friction pressure values do not match at all. The static friction slip vectors 
were examined, and they are not perfectly vertical, and they are not fully 
loaded. This results in the static friction force being less than the maximum 
possible, and thus the pressure is higher since less of the force is removed by 
the walls. 
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Figure 5.8: Comparison With 0.05 and 0.15 /i 
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Figure 5.9: Comparison With 0.25 and 0.9 /i 



Chapter 6 

Unphysical Approximations 



Modeling the full physical effects that occur in a pebble bed reactor mechan- 
ics is not computationally possible with current computer resources. In fact, 
even modeling all the intermolecular forces that occur between two pebbles 
at sufficient levels to reproduce all macroscopic behavior is probably com- 
putationally intractable at the present time. This is partially caused by the 
complexity of effects such as inter-grain boundaries and small quantities of 
impurities that affect the physics and different levels between the atomic 
effects and the macroscopic world. Instead, all attempts at modeling the 
behavior of pebble bed reactor mechanics have relied on approximation to 
make the task computationally practical. The PEBBLES simulation has as 
high or higher fidelity than past efforts, but it does use multiple unphysical 
approximations. This chapter will discuss the approximations so that future 
simulation work can be improved, and an understanding of what limitations 
exist when applying PEBBLES to different problems. 

In different regions of the reactor, the radioactivity and the fission will 
heat the pebbles differently, and the fiow of the coolant helium will distribute 
this heat around the reactor. This will change the temperature of different 
parts of the reactor. Since the temperature will be different, the parameters 
driving the mechanics of the pebbles will be different as well. This includes 
parameters such as the static friction coefficients and the size of the pebbles 
which will change through thermal expansion. As well, parameters such as 
static friction can also vary depending on the gas in which they currently are 
in and in which they were, since some of the gas tends to remain in and on 
the carbon surface. Graphite dust produced by wear may also affect static 
friction in downstream portions of the reactor. 
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The pebbles in a pebble bed reactor have helium gas flowing around and 
past them. PEBBLES and all other pebble bed simulations ignore effects 
of this on pebble movement. However, the gas will cause both additional 
friction when the pebbles are dropping through the reactor, and the motion 
of gas will cause additional forces on pebbles. 

Pebble bed mechanics simulations use soft spheres. Physically, there will 
be deflection of spheres under pressure (even the pressure of just one sphere 
on the floor), but the true compression is much smaller than what is actu- 
ally modeled. In PEBBLES, the forces are chosen to keep the compression 
distance at a millimeter or below. Another effect related to the physics sim- 
ulation is that force is transmitted via contact. This means the force from 
one end of the reactor is transmitted at a speed related to the time-step used 
for the simulation, instead of the speed of sound. 

Since simulating billions of time-steps is time consuming, two approxi- 
mations are made. First, instead of simulating the physical time that pebble 
bed reactors have between pebble additions (on the order of 2-5 minutes), 
new pebbles are added at a rate between a quarter second and two seconds. 
This may result in somewhat unphysical simulations since some vibration 
that would have dampened out with a longer time between pebble additions 
still exists when the next pebble impacts the bed. Second, since full recircu- 
lation of all the pebbles is computationally costly, for some simulations, only 
a partial recirculation or no recirculation is done. 

The physics models do not take into account several physical phenomena. 
The physics do not handle pure spin effects, such as when two pebbles are 
contacting and are spinning with an axis around the contact point. This 
should result in forces on the pebbles, but the physics model does not handle 
this effect since the contact velocity is calculated as zero. In addition, when 
the pebble is rolling so that the contact velocity is zero because the pebble's 
turning axis is parallel to the surface and at the same rate as the pebble is 
moving along the surface, there should be rolhng friction, but this effect is 
not modeled either. As well, the equations used assume that the pebbles are 
spherically symmetric, but defects in manufacturing and slight asymmetries 
in the TRISO particle distribution mean that there will be small deviations 
from being truly spherically symmetric. 

The physics model does not match classical Hertzian or Mindlin and Dere- 
siewicz elastic sphere behavior. The static friction model is a simpliflcation 
and does not capture all the hysteretic effects of true static friction. Effec- 
tively, this means that hg, the coefficient used to calculate the force from 
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slip, is not a constant. In order to fully discuss this, some features of these 
models will be discussed in the following paragraphs. 

Since closed-form expressions exist for elastic contact between spheres, 
they will be used, instead of a more general case which lacks closed-form 
expressions. Spheres are not a perfect representation of the effect of contact 
between shapes such as a cone and a sphere, but should give an approxima- 
tion of the size of the effect of curvature. 

The amount of contact area and displacement of distant points for two 
spheres or one sphere and one spherical hole (that is negative curvature) for 
elastic spheres can be calculated via Hertzian theory ( Johnson| |1985 ). For 
two spherical surfaces the following variables are defined: 

R Ri R2 ^ ^ 

and 

1 1 - z/2 1 - zy2 

— = + 6.2 

E* El E2 ^ ' 

with Ri the ith's sphere's radius, Ei the Young's modulus, z/, the Poisson's 
ration of the material. For a concave sphere, the radius will be negative. 
Then, via Hertzian theory, the contact circle radius will be: 
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(6.3) 



where is the normal force. The mutual approach of distant points is 
given by: 
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Notice that the above differs compared to the Hooke's Law formulation 
that PEBBLES uses. The maximum pressure will be: 
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So as a function of the radii -Ri and R2, the circle radius of the contact 
will be: 
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(6.6) 



If m is used for the multiple of negative curvature sphere of the radius of 
the other, then the equation becomes: 
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which can be rearranged to: 
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From this equation, as m increases, it has less effect on contact area, so 
if R2 is much greater than Ri, the contact area will tend to be dominated 
by Ri rather than R2. For example, typical radii in PEBBLES might be 18 
cm outlet chute and a 3 cm pebble, which would put m at 6, so the effect 
on contact area radius would be about 33% difference compared to pebble 
to pebble contact area radius, or 6% compared to a fiat surface]^ 

To some extent, the actual contact area is irrelevant for calculating the 
maximum static friction force as long as some conditions are met. Both 
surfaces need to be of a uniform material. The basic macroscopic description 
IF5I <= yu|A^| needs to hold, so changing the area changes the pressure 
P = N/a, but not the maximum static friction force. If the smaller area 
causes the pressure to increase enough to cause plastic rather than elastic 
contact, then through that mechanism, the contact area would cause actual 
differences in experimental values. PEBBLES also does not calculate plastic 
contact effects. 

The contact area causes an effect through another mechanism. The tan- 
gential compliance in the case of constant normal and increasing tangential 



1 



-1/3. 



^ Sample values of A: = (l — ^ 
6, fc = 0.94 sphere to outlet chute and m 



m = — l.fc = 1.26 for sphere to sphere, m 
= 00, /c = 1 sphere to flat plane. 
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force, that is the slope of the curve relating displacement to tangential force, 
is given in Mindlin and Deresiewicz as: 

2 -i^ . . 

6.9 

Sua ^ ' 

Since the contact area radius, a, is a function of curvature, the slope of the 
tangential compliance will be as well, which is another effect that PEBBLES' 
constant hs does not capture. 

In summary for the static friction using a constant coefficient for hs yields 
two different approximations. First, using the same constants for wall con- 
tact when there is different curvatures is an approximation that will give 
somewhat inconsistent results. Since the equations for spherical contact are 
dominated by the smaller radius object, this effect is somewhat less but still 
exists. Second, using the same constant coefficient for different loading his- 
tories is a approximation. For a higher fidelity, these effects need to be taken 
into account. 



Chapter 7 



Code Speedup and 
Parallelization 

Planned and existing pebble bed reactors can have on the order of 100,000 
pebbles. For some simulations, these pebbles need to be followed for long 
time periods, which can require computing billions of time-steps. Multiply- 
ing the time-steps required by the number of pebbles being computed over 
leads to the conclusion that large numbers of computations are required. 
These computations should be as fast as possible, and should be as parallel 
as possible, so as to allow relevant calculations to be done in a reasonable 
amount of time. This chapter discusses the process of speeding up the code 
and parallelizing it. 

The PEBBLES program has three major portions of calculation. The first 
is determining which pebbles are in contact with other pebbles. The second 
computational part is determining the time derivatives for all the vectors for 
all the pebbles. The third computational part is using the derivatives to up- 
date the values. Overall, for calculation of a single time-step, the algorithm's 
computation time is linearly proportional to the number of pebbles, that is 
0(nH 



^O(n): The algorithm scales linearly (n) with increasing input size. So if it runs with 
100 pebbles it takes roughly 10 times as long as when it runs it only 10 pebbles. Or if it 
goes from 10 pebbles to 20 pebbles it will take twice as long to run. 
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7.1 General Information about profiling 

There are four different generic parts of the complete calculation that need 
to considered for determining the overall speed. The first consideration is 
the time to compute arithmetic operations. Modern processors can com- 
plete arithmetic operations in nanoseconds or fractions of nanoseconds. In 
the PEBBLES code, the amount of time spent on arithmetic is practically 
undetectable in wall clock changes. The second consideration is the time re- 
quired for reading memory and writing memory. For main memory accesses, 
this takes hundreds of CPU clock cycles, so these times are on the order of 



fractions of microseconds (Drepper 2007). Because of the time required to 



access main memory, all modern CPUs have on-chip caches, that contain a 
copy of the recently used data. If the memory access is in the CPU's cache, 
the data can be retrieved and written in a small number of CPU cycles. Main 
memory writes are somewhat more expensive than main memory reads, since 
any copies of the memory that exist in other processor's caches need to be 
updated or invalidated. So for a typical calculation like a + 6 — )■ c the time 
spent doing the arithmetic is trivial compared to the time spent reading in 
a and b and writing out c. 

The third consideration is the amount of time required for parallel pro- 
gramming constructs. Various parallel synchronization tools such as atomic 
operations, locks and critical sections take time. These take an amount of 
time on the same order of magnitude as memory writes. However, they typi- 
cally need a read and then a write without any other processor being able to 
access that chunk of memory in between which requires additional overhead, 
and a possible wait if the memory address is being used by another process. 
Atomic operations on x86_64 architectures are faster than using locks, and 
locks are generally faster than using critical sections. The fourth consider- 
ation is network time. Sending and receiving a value can easily take over a 
millisecond for the round trip time. These four time consuming operations 
need to be considered when choosing algorithms and methods of calculation. 

There are a variety of methods for profiling the computer code. The 
simplest method is to use the FORTRAN 95 intrinsics CPU_TIME and DATE_- 
AND_TIME. The CPU_TIME subroutine returns a real number of seconds of CPU 
time. The DATE_AND_TIME subroutine returns the current wall clock time in 
the VALUES argument. With gfortran both these times are accurate to 
at least a millisecond. The difference between two different calls of these 
functions provide information on both the wall clock time and the CPU 
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time between the calls. (For the DATE_AND_TIME subroutine, it is easiest 
if the days, hours, minutes, seconds and milliseconds are converted to a 
real seconds past some arbitrary time.) The time methods provide basic 
information and a good starting point for determining which parts of the 



program are consuming time. For more detailed profiling the oprof ile (opr 



2009 ) program can be used on Linux. This program can provide data at the 



assembly language level which makes it possible to determine which part of 
a complex function is consuming the time. Non-assembly language profilers 
are difficult to accurately use on optimized code, and profiling non-optimized 
code is misrepresentative. 



7.2 Overview of Parallel Architectures and 
Coding 

Parallel computers can be arranged in a variety of ways. Because of the 
expense of linking shared memory to all processors, a common architecture 
is a cluster of nodes with each node having multiple processors. Each node 
is linked to other nodes via a fast network connection. The processors on 



a single node share memory. Figure 7.1 shows this arrangement. For this 
arrangement, the code can use both the OpenMP (Open Multi-Processing) 
(ope, 2008) and the MPI (Message Passing Interface) ( |mpi 2009) libraries. 
MPI is a programming interface for transferring data across a network to 
other nodes. OpenMP is a shared memory programming interface. By us- 
ing both programming interfaces high speed shared memory accesses can be 
used on memory shared on the node and the code can be parallelized across 
multiple nodes. 
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Figure 7.1: Sample Cluster Architecture 



42 



CHAPTER 7. CODE SPEEDUP AND PARALLELIZATION 



7.3 Lock- less Parallel 0(N) Collision Detec- 
tion 

For any granular material simulation, which particles are in contact must be 
determined quickly and accurately for each time-step. This is called collision 
detection, though for pebble simulations it might be more accurately labeled 
contact detection. The simplest algorithm for collision detection is to iterate 
over all the other objects and compare each one to the current object for 
collision. To determine all the collisions using that method, 0{N'^) time is 
required. 



An improved algorithm by Cohen et al. (1995) uses six sorted lists of 



the lower and upper bounds for each object. (There is one upper bound 
list and one lower bound list for each dimension.) With this algorithm, to 
determine the collisions for a given object, the bounds of the current objects 
are compared to bounds in the list — only objects that overlap the bounds 
in all three dimensions will potentially collide. This algorithm typically has 
approximately 0{N \og[N)) timen because of the sorting of the bounding 



lists (Cohen et al. 1995) 



A third, faster method, grid collision detection, is available if the following 
requirements hold: (1) there is a maximum diameter of object, and no object 
exceeds this diameter, and (2) for a given volume, there is a reasonably 
small, finite, maximum number of objects that could ever be in that volume. 
These two constraints are easily satisfied by pebble bed simulations, since the 
pebbles are effectively the same size (small changes in diameter occur due to 
wear and thermal effects). A three-dimensional parallelepiped grid is used 
over the entire range in which the pebbles are simulated. The grid spacing 
gs is set at the maximum diameter of any object (twice the maximum radius 
for spheres). 

Two key variables are initialized, grid-Count{x,y, z), the number of peb- 
bles in grid locations x,y,z; and gridJds{x,y, z,i), the pebble identification 
numbers (ids) for each x,y,z location. The id is a unique number assigned to 
each pebble in the simulation. The spacing between successive grid indexes 
is gs, so the index of a given x location can be determined by [{x — Xmin)/ gs\ 
where Xmin is the zero x index's floor; similar formulas are used for y and z. 

The grid is initialized by setting grid.count{:, :, :) = 0, and then the x,y,z 
indexes are determined for each pebble. The grid-Count at that location is 



^In order from slowest to fastest (for sufficiently big N): 0{N'^),0{NlogiN),0{N),0{l). 



7.3. LOCK-LESS PARALLEL 0(N) COLLISION DETECTION 



43 



then atomicalljjj incremented by one and fetched. Because OpenMP 3.0 does 
not have a atomic add-and-fetch, the lock xadd x86_64 assembly language 
instruction is put in a function. The grid.count provides the fourth index 
into the gridJds array, so the pebble id can be stored into the ids array. 
The amount of time to zero the grid_count array is a function of the volume 
of space, which is proportional to the number of pebbles. The initialization 
iteration over the pebbles can be done in parallel because of the use of an 
atomic add-and-fetch function. Updating the grid iterates over the entire list 
of pebbles so the full algorithm for updating the grid is 0{N) for the number 
of pebbles. 

Once the grid is updated, the nearby pebbles can be quickly determined. 
Figure [7!2l illustrates the general process. First, index values are computed 
from the pebble and used to generate xc, yc, and zc. This finds the center 
grid location, which is shown as the bold box in the figure. Then, all the 
possible pebble collisions must have grid locations (that is, their centers are 
in the grid locations) in the dashed box, which can be found by iterating 
over the grid locations from xc — 1 to xc + 1 and repeating for the other two 
dimensions. There are 3^ grid locations to check, and the number of pebbles 
in them are bounded (maximum 8), so the time to do this is bounded. Since 
this search does not change any grid values, it can be done in parallel without 
any locks. 




XC 



Figure 7.2: Determining Nearby Pebbles from Grid 

Therefore, because of the unique features of pebble bed pebbles simu- 
lation, a parallel lock- less 0{N) algorithm for determining the pebbles in 
contact can be created. 



■^In this chapter, atomic means uncutable, that is the entire operation is done in one 
action without interference from other processors. 
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7.4 MPI Speedup 

The PEBBLES code uses MPI to distribute tiie computational work across 
different nodes. Tfie MPI/OpenMP fiybrid parallelization splits the calcu- 
lation of the derivatives and the new variables geometrically and passes the 
data at the geometry boundaries between nodes using messages. Each peb- 
ble has a primary node and may also have various boundary nodes. The 
pebble-primary-node is responsible for updating the pebble position, veloc- 
ity, angular velocity, and slips. The pebble-primary-node also sends data 
about the pebble to any nodes that are the pebble boundary nodes and will 
transfer the pebble to a different node if the pebble crosses the geometric 
boundary of the node. Boundary pebbles arc those close enough to a bound- 
ary that their data needs to be present in multiple nodes so that the node's 
primary pebbles can be properly updated. Node is the master node and 
does processing that is simplest to do on one node, such as writing restart 
data to disk and initializing the pebble data. The following steps are used 
for initializing the nodes and then transferring data between them: 

1. Node calculates or loads initial positions of pebbles. 

2. Node creates the initial domain to node mapping. 

3. Node sends domain to node mapping to other nodes. 

4. Node sends other nodes their needed pebble data. 

Order of calculation and data transfers in main loop: 

1. Calculate derivatives for node primary and boundary pebbles. 

2. Apply derivatives to node primary pebble data. 

3. For every primary pebble, check with the domain module to determine 
the current primary node and any boundary nodes. 

(a) If the pebble now has a different primary node, add the pebble id 
to the transfer list to send to the new primary node. 

(b) If the pebble has any boundary nodes, add the pebble id to the 
boundary send list to send it to the node for which it is a boundary. 
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4. If this is a time step where Node needs all the pebble data (such as 
when restart data is being written), add all the primary pebbles to the 
Node boundary send list. 

5. Send the number of transfers and the number of boundary sends that 
this node has to all the other nodes using buffered sends. 

6. Initialize three Boolean lists of other nodes that this node has: 

(a) data-to-send-to with "true" if the number of transfers or boundary 
sends is nonzero, and "false" otherwise 

(b) received-data-from to "false" 

(c) received-the-number-of-transfers and the number-of-boundary-sends 
with "false." 

7. While this node has data to send to other nodes and other nodes have 
data to send to this node loop: 

(a) Probe to see if any nodes that this node needs data from have 
data available. 

i. If yes, then receive the data and update pebble data and the 
Boolean lists as appropriate 

(b) If there are any nodes that this node has data to send to, and this 
node has received the number of transfers and boimdary sends 
from, then send the data to those nodes and update the Boolean 
data send hst for those nodes. 

8. Flush the network buffers so any remaining data gets received. 

9. Node calculates needed tallies. 

10. If this is a time to rebalance the execution load: 

(a) Send wall clock time spent computing since last rebalancing to 
node 

(b) Node uses information to adjust geometric boundaries to move 
work towards nodes with low computation time and away from 
nodes with high computation time 
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(c) Node sends new boundary information to other nodes, and 
needed data to other nodes. 

11. Continue to next time step and repeat this process until all time-steps 
have been run. 

All the information and subroutines needed to calculate the primary and 
boundary nodes that a pebble belongs to are calculated and stored in a FOR- 
TRAN 95 module named network^domainjnodule. The module uses two 
derived types: network_domain_type and network_domain_location_type. 
Both types have no public components so the implementation of the domain 
calculation and the location information can be changed without changing 
anything but the module, and the internals of the module can be changed 
without changing the rest of the PEBBLES code. The location type stores 
the primary node and the boundary nodes of a pebble. The module contains 
subroutines for determining the location type of a pebble based on its posi- 
tion, primary and boundary nodes for a location type, and subroutines for 
initialization, load balancing, and transferring of domain information over 
the network. 

The current method of dividing the nodes into geometric domains uses a 
list of boundaries between the z (axial) locations. This list is searched via 
binary search to find the nodes nearest to the pebble position, as well as 
those within the boundary layer distance above and below the zone interface 
in order to identify all the boundary nodes that participate in data transfers. 
The location type resulting from this is cached on a fine grid, and the cached 
value is returned when the location type data is needed. The module contains 
a subroutine that takes a work parameter (typically, the computation time 
of each of the nodes) and can redistribute the z boundaries up or down to 
shift work towards nodes that are taking less time computing their share 
of information. If needed in the future, the z-only method of dividing the 
geometry could be replaced with a full 3-D version by modifying the network 
domain module. 

7.5 OpenMP Speedup 

The PEBBLES code uses OpenMP to distribute the calculation over multiple 
processes on a single node. OpenMP allows directives to be given to the 
compiler that direct how portions of code are to be parallelized. This allows 
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a single piece of code to be used for both the single processor version and 
the OpenMP version. The PEBBLES parallelization typically uses OpenMP 
directives to cause loops that iterate over all the pebbles to be run in parallel. 

Some details need to be taken into consideration for the parallehzation 
of the calculation of acceleration and torque. The physical accelerations 
imposed by the wall are treated in parallel, and there is no problem with 
writing over the data because each processor is assigned a portion of the 
total zone inventory of pebbles. For calculating the pebble-to-pebble forces, 
each processor is assigned a fraction of the pebbles, but there is a possibility 
of the force addition computation overwriting another calculation because 
the forces on a pair of pebbles are calculated and then the calculated force 
is added to the force on each pebble. In this case, it is possible for one 
processor to read the current force from memory and add the new force 
from the pebble pair while another processor is reading the current force 
from memory and adding its new force to that value; they could both then 
write back the values they have computed. This would be incorrect because 
each calculation has only added one of the new pebble pair forces. Instead, 
PEBBLES uses an OpenMP ATOMIC directive to force the addition to be 
performed atomically, thereby guaranteeing that the addition uses the latest 
value of the force sum and saves it before a different processor has a chance 
to read it. 

For calculating the sum of the derivatives using Euler's method, updating 
concurrently poses no problem because each individual pebble has derivatives 
calculated. The data structure for storing the pebble-to-pebble slips (sums 
of forces used to calculate static friction) is similar to the data structure used 
for the coUision detection grid. A 2-D array exists where one index is the 
from-pebble and the other index is for storing ids of the pebbles that have slip 
with the first pebble. A second array exists that contains the number of ids 
stored, and that number is always added and fetched atomically, which allows 
the slip data to be updated by multiple processors at once. These combine 
to allow the program to run efficiently on shared memory architectures. 

7.6 Checking the Parallelization 

The parallelization of the algorithm is checked by running the test case with 
a short number of time steps (10 to 100). Various summary data are checked 
to make sure that they match the values computed with the single processor 



48 



CHAPTER 7. CODE SPEEDUP AND PARALLELIZATION 



version and between different numbers of nodes and processors. For example, 
with the NGNP-600 model used in the results section, the average overlap of 
pebbles at the start of the run is 9.665281e-5 meters. The single processor 
average overlap at the end of the 100 time-step run is 9.693057e-5 meters, 
the 2 nodes average overlap is 9.693043e-5 meters, and the 12 node average 
overlap is 9.693029e-5 meters. The lower order numbers change from run to 
run. The start-of-run values match each other exactly, and the end-of-run 
values match the start of run values to two significant figures. However, the 
three different end-of-run values match to five significant digits. In short, the 
end values match each other more than they match the start values. The 
overlap is very sensitive to small changes in the calculation because it is a 
function of the difference between two positions. During coding, multiple de- 
fects were found and corrected by checking that the overlaps matched closely 
enough between the single processor calculation and the multiple processor 
calculations. The total energy or the linear energy or other computations 
can be used similarly since the lower significant digits also change frequently 
and are computed over all the pebbles. 



7.7 Results 



The data in Table 7_A and Table 7.2| provide information on the time used 



with the current version of PEBBLES for running 80 simulation time steps 
on two models. The NGNP-600 model has 480,000 pebbles. The AVR model 
contains 100,000 pebbles. All times are reported in units of wall-clock sec- 
onds. The single processor NGNP-600 model took 251 seconds and the AVR 
single processor model took 48 seconds when running the current version. 
The timing runs were carried out on a cluster with two Intel Xeon X5355 
2.66 GHz processors per node with a DDR 4X InfiniBand interconnect net- 
work. The nodes had 8 processors per node. The gfortran 4.3 compiler was 
used. 

Significant speedups have resulted with both the OpenMP and MPI- 
/OpenMP versions. A basic time step for the NGNP-600 model went from 
3.138 seconds to 146 milliseconds when running on 64 processors. Since a full 
recirculation would take on the order of 1.6e9 time steps, the wall clock time 
for running a full recirculation simulation has gone from about 160 years to 
a little over 7 years. For smaller simulation tasks, such as simulating the mo- 
tion of the pebbles in a pebble bed reactor during an earthquake, the times 
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Processes 



Table 7.1: O 
AVR Speedup 



aenMP speedup results 
Efficiency NGNP-600 Speedup 



Efficiency 



47.884 
53.422 
29.527 
21.312 
16.660 
13.884 
12.012 
10.698 
9.530 



1 

0.89633 
1.6217 
2.2468 
2.8742 
3.4489 
3.98635 
4.4760 
5.0246 



100.00% 

89.63% 

81.09% 

74.89% 

71.85% 

68.98% 

66.44% 

63.94% 

62.81% 



251.054 

276.035 

152.479 

104.119 

80.375 

68.609 

61.168 

54.011 

49.171 



1 

0.90950 
1.6465 
2.4112 
3.1235 
3.6592 
4.1043 
4.6482 
5.1057 



100.00% 

90.95% 

82.32% 

80.37% 

78.09% 

73.18% 

68.41% 

66.40% 

63.82% 



Table 7.2: MPI/OpenMP speedup results 



Nodes 


Procs 


AVR 


Speedup 


Efficiency 


NGNP-600 


Speedup 


Efficiency 


1 


1 


47.884 


1 


100.00% 


251.054 


1 


100.00% 


1 


8 


10.696 


4.4768 


55.96% 


55.723 


4.5054 


56.32% 


2 


16 


6.202 


7.7207 


48.25% 


30.642 


8.1931 


51.21% 


3 


24 


4.874 


9.8244 


40.93% 


23.362 


10.746 


44.78% 


4 


32 


3.935 


12.169 


38.03% 


17.841 


14.072 


43.97% 


5 


40 


3.746 


12.783 


31.96% 


16.653 


15.076 


37.69% 


6 


48 


3.534 


13.550 


28.23% 


15.928 


15.762 


32.84% 


7 


56 


3.285 


14.577 


26.03% 


15.430 


16.271 


29.05% 


8 


64 


2.743 


17.457 


27.28% 


11.688 


21.480 


33.56% 


9 


72 


2.669 


17.941 


24.92% 


11.570 


21.699 


30.14% 


10 


80 


2.657 


18.022 


22.53% 


11.322 


22.174 


27.72% 


11 


88 


2.597 


18.438 


20.95% 


11.029 


22.763 


25.87% 


12 


96 


2.660 


18.002 


18.75% 


11.537 


21.761 


22.67% 



are more reasonable, taking about 5e5 time steps. Thus, for the NGNP-600 
model, a full earthquake can be simulated in about 20 hours when using 64 
processors. For the smaller AVR model, the basic time step takes about 34 
milliseconds when using 64 processors. Since there are less pebbles to recir- 
culate, a full recirculation would take on the order of 2.5e8 time steps, or 
about 98 days of wall clock time. 



Chapter 8 
Applications 



The knowledge of the packing and flow patterns (and to a much lesser extent 
the position) of pebbles in the pebble bed reactor is an essential prerequisite 
for many in-core fuel cycle design activities as well as for safety assessment 
studies. Three applications have been done with the PEBBLES code. The 
major application is the computation of pebble positions during a simulated 
earthquake. Two other applications that have been done are calculation of 
space dependent Dancoff factors and calculation of the angle of repose for a 
HTR-10 simulation. 



8.1 Applications in Support of Reactor Physics 
8.1.1 DancofF Factors 

The calculation of Dancoff factors is an example application that needs ac- 
curate pebble position data. The Dancoff factor is used for adjusting the 
resonance escape probability for neutrons. There are two Dancoff factors 
that use pebble position data. The first is the inter-pebble Dancoff factor 
that is the probability that a neutron escaping from the fuel zone of a pebble 
crosses a fuel particle in another pebble. The second is the pebble-pebble 
Dancoff factor, which is the probability that a neutron escaping one fuel zone 
will enter another fuel zone without interacting with a moderator nuclide. 



Kloosterman and Ougouag (2005) use pebble location information to calcu- 
late the probability by ray tracing from fuel lumps until another is hit or the 
ray escapes the reactor. The PEBBLES code has been used for providing 
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average velocity 



/ / 
/ / 



0.1 0,2 0.3 



Figure 8.1: Flow Field Representation (arrow lengths are propor- 
tional to local average pebble velocity) 



position information to J. L. Kloosterman and A. M. Ougouag's PEBDAN 



program. This program calculate these factors as shown in Figure 8.2 which 
calculates them for the AVR reactor model. 



8.1.2 Angle of Repose 

The PEBBLES code was used for calculating the angle of repose for an 



analysis of the HTR-10 first criticality (Terry et al. 2006). The pebble bed 



code recirculated pebbles to determine the angle at which the pebbles would 



stack at the top of the reactor as shown in Figure 8.3, since this information 



was not provided, but was needed for the simulation of the reactor (Baumer 



et al. 1990). 



8.1.3 Pebble Ordering with Recirculation 

During experimental work before the construction of the AVR, it was dis- 
covered that when the pebbles were recirculated, the ordering in the pebbles 
increased. Figures 8.4 and 8.5 show that this effect occurs in the PEBBLES 
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Dancoff Factors 




Figure 8.2: Dancoff Factors for AVR 

simulation as well. The final AVR design incorporated indentations in the 
wall to prevent this from occurring. 

8.2 Application to Earthquake modeling 

The packing fraction of the pebbles in a pebble bed reactor can vary depend- 
ing on the method of packing and the subsequent history of the packing. 
This packing fraction can affect the neutronics behavior of the reactor, since 
it translates into an effective fuel density. During normal operation, the 
packing fraction will vary only slowly, over the course of weeks and then 
stabilize. During an earthquake, the packing fraction can increase suddenly. 
This packing fraction change is a concern since packing fraction increase can 
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Figure 8.3: 2-D Projection of Pebble Cone on HTR-10 (crosses rep- 
resent centers of pebbles) 



Figure 8.4: Pebbles Before Recirculation 



increase the neutron multiplication and cause criticality concerns as shown 
by Ougouag and Terry (2001). 

The PEBBLES code can simulate this increase and determine the rate of 
change and the expected final packing fraction, thus allowing the effect of an 
earthquake to be simulated. 
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Figure 8.5: Pebbles After Recirculation 



8.2.1 Movement of Earthquakes 



The movement of earthquakes has been well studied in the past. The mag- 
nitude of the motion of earthquakes is described by the Mercalli scale, which 
describes the maximum acceleration that a given earthquake will impart to 
structures. For a Mercalli X earthquake, the maximum acceleration is about 
1 g. The more familiar Richter scale measures the total energy release of an 
earthquake (Lamarsh, 1983), which is not useful for determining the effect on 
a pebble bed core. For a given location, the soil properties can be measured, 
and using soil data and the motion that the bedrock will undergo, the motion 
on the surface can be simulated. The INL site had this information gener- 
ated in order to determine the motion from the worst earthquake that could 
be expected over a 10,000 years period (Payne, 2003). This earthquake has 
roughly a Mercalli IX intensity. The data for such a 10,000 year earthquake 
are used for the simulation in this dissertation. 
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8.2.2 Method Of Simulation 

The code simulates earthquakes by adding a displacement to the walls of the 
reactor. As well, the velocity of the walls needs to be calculated. The dis- 
placement in the simulation can be specified either as the sum of sine waves, 
or as a table of displacements that specifies the x, y, and z displacements for 
each time. At each time step both the displacement and the velocity of the 
displacement are calculated. When the displacement is calculated by a sum 
of sine functions, the current displacement is calculated by adding vector di- 
rection for each wave and the velocity is calculated from the sum of the first 
derivative of all the waves. When the displacement is calculated from a table 
of data, the current displacement is a linear interpolation of the two nearest 
data points in the table, and the velocity is the slope between them. The 
walls are then assigned the appropriate computed displacement and velocity. 
Figure |8.6| shows the total displacement for the INL earthquake simulation 
specifications that were used in this paper. 



8.2.3 Earthquake Results 

The results of two simulations carried out here show a substantially safer 



behavior than the Ougouag and Terry (2001) bounding calculations. The 



methodology was applied to a model of the PBMR-400 model and two dif- 
ferent static friction coefficients were tested, 0.65 and 0.35. The packing 
fraction increased from 0.593 to 0.594 over the course of the earthquake with 
the 0.65 static friction model, with the fastest increase was from 0.59307 to 
0.59356 and took place over 0.8 seconds. With the 0.35 static friction model, 
the overall increase was from 0.599 to 0.601. The fasted increase was from 
0.59964 to 0.60011 in 0.8 seconds. This is remarkably small when compared 
to the bounding calculation packing fraction increase rate of 0.129 sec~^ in 
free fallj^ Both computed increases and packing fraction change rates are 
substantially below the free fall bounding rate and packing fraction change 
of a transition from 0.60 to 0.64 in 0.31 seconds. The computed rate and the 
total packing fraction increase are in the range that can be offset by thermal 
feedback effects for uranium fueled reactors. 

During the course of the earthquake, the boundary density fluctuations 
(that is the oscillations in packing fraction near a boundary) are observed 



^The free fall rate is determined by calculating the packing fraction increase if the 
pebbles were in gravitational free fall. 
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Figure 8.6: Total Earthquake Displacement 



to increase in amplitude. Figure 8.9 shows the packing fraction before the 
earthquake and after the earthquake in the radial direction. These were taken 
from 4 to 8 meters above the fuel outlet chute in the PBMR-400 model. 
All the radial locations have increased packing compared to the packing 
fraction before the earthquake, but the points that are at boundary density 



fluctuation peaks increase the most. This effect can be seen in figure 8.10 



which shows the increase in packing fraction before the earthquake and after 

A previous version of the positional data from the earthquake simulation 
was provided to J. Ortensi. This data was used by him to simulate the effects 
of an earthquake on a pebble bed reactor (Ortensi, 2009). Essentially, two 



factors cause an increase in reactivity. The first is the increased density of 
the pebbles and the second is due to the control rods being at the top of 
reactor, so when the top pebbles move down the control rod worth (effect) 



8.2. APPLICATION TO EARTHQUAKE MODELING 



57 



o 



o 
CO 



0.596 



0.595 



0.594 



0.593 - 



-g 0.592 
as 

a. 



0.591 



0.59 



0.589 





2to3 
3to4 
4to5 
5to6 
6to7 
7to8 
8to9 
9tol0 
lOtoll 
lltol2 
I 



-B- 



-e- 



10 15 20 25 
Time (seconds) 



30 35 40 45 



Figure 8.7: 0.65 Static Friction Packing over Time 



decreases. However, the reactivity increase causes the fuel temperature to 
rise, which causes Doppler broadening and more neutrons are absorbed by 



the ^^^U, which causes the reactivity to fall. Figure 8.11 shows an example 
of this. 



8.2.4 Earthquake Equations 

For each time-step, the simulation calculates both a displacement and a wall 
velocity. 

For the sum of waves method, the displacement is calculated by: 



Ed 



sin ( [t — bj h c I + o 



•1) 
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Figure 8.8: 0.35 Static Friction Packing over Time 

where t is the current time, S is the time the wave starts, p is the period 
of the wave, c is the initial cycle of the wave, o is the offset, and D is the 
maximum displacement vector. 

The velocity is calculated by: 



m 



27rD A ^,27r 
5^— CO. {^{t-S)j + c 



(8.2) 



For the tabular data, the displacement and velocity are calculated by: 



d = (1 - o)Tk + oTfc+i 
m = \{Tk+i - Tfc) 





(8.3) 
(8.4) 
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Figure 8.9: Different Radial Packing Fractions 

where Tj is the displacement at the ith time-step, o is a number between 
and 1 that specifies where is the start of the time-step and 1 is the end, 
and 6 is the time in seconds between time-steps. 

With these displacements, the code then uses: 

p = p + d (8.5) 
v' = V + m (8.6) 

as the adjusted position and velocity. 
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Figure 8.10: Changes in Packing Fraction 
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Figure 8.11: Neutronics and Thermal Effects from J. Ortensi 



Chapter 9 



Construction of a Dust 
Production Framework 



With the creation of the PEBBLES simulation, one issue that was examined 
was using the simulation to attempt to predict the volume of dust that would 
be produced by an operating pebble bed reactor. This is an important issue 
that could affect the choice of a pebble bed reactor versus a prismatic reactor 
for process heat applications. However, as this chapter and Appendix [B] will 
discuss, while the PEBBLES code has the force and motion data required 
for this simulation, the coefficients that would allow this information to be 
used have not been sufficiently robustly experimentally determined yet. 

With the data provided by PEBBLES, equations to link the dust pro- 
duction to PEBBLES calculated quantities were examined. As shown in 
equation B.l the volume of dust produced can be approximated if the nor- 
mal force of contact, the length slide and the wear coefficients are known. 
The force of contact and the length slide are calculated as part of the PEB- 
BLES simulation, so this method was used to calculate dust production for 
the AVR reactor. This resulted in an estimate of four grams of graphite dust 
produced per year as compared to the measured value of three kilograms of 
dust produced per year. Several possible causes of this were identified in the 
paper documenting this work (Cogliati and Ougouag, 2008). A key first is- 
sue as described by this dissertation's previous work section is that there are 
no good measurements of graphite wear rates in pebble bed reactor relevant 
conditions (especially for a reactor run at AVR temperatures). A second 
issue is that the previous model of AVR was missing features including the 
reactor control rod nose cones and wall indentations. A third issue, identi- 
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Figure 9.1: Pebble to Pebble Distance Traveled 



fied after the paper's publication, is that significant portions of the length 
traveled were due not to motion down through the reactor. Instead, much 
of the length that was tallied was due to pebbles vibrating. In the model 
used in the paper, this problem was about four times more severe than the 



current model, due to the new addition of slip correction via equation 5.2 



As an illustration of the general framework for dust production, a simple 
cylindrical vat the size of AVR was simulated. In this model, the outlet chute 
starts shrinking at height zero. The length (L), and length times normal force 



(NL) were tallied on intervals of 6 cm and the figures 9.1 to 9.5 show the 



results after recirculating 400 pebbles. The results are on a per pebble-pass 
basis. Two sets of static friction and kinetic friction pairs are used, one with 
a static friction coefficient of 0.35 and kinetic of 0.25, and the other with 
static of 0.65 and kinetic of 0.4. Figure 9.1 shows the calculated pebble-to- 
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Figure 9.2: Pebble to Surface Distance Traveled 



pebble lengths. Notice that for the 0.65 static friction simulation, about 10 
meters of length traveled is occurring at the peak in a 6 cm long tally. Since 
the pebble is traveling about 0.06 m and has at most 12 pebble to pebble 
contacts, essentially all but a small portion of this length traveled is due to 
when the pebbles vibrate relative to each other. This vibration is caused by 
the impact of the pebbles coming from the inlet chutes and hitting the top 
of the bed. This likely is a true physical effect, which has not been discussed 
in literature this author is aware of. However, in order to obtain the correct 
magnitude of this vibrational effect, two things must be correct. First, the 
simulation must dissipate the vibration at the correct rate, and second, the 
pebbles must impact the bed at the correct velocity. Quality estimates of 
both will need to be made to finish the dust production work. Note that 
with the lower static friction value, the effect is even more pronounced. 
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Figure 9.3: Average Normal Contact Force 



Figure 9.2 is also expected to have a vibrational component, since only 
a small portion of pebbles should contact the wall, and therefore for a 6 cm 
tally, the length traveled by the average pebble should be much lower than 6 
cm. As the chute is entered, the distance the average pebble travels increases 
in the 0.65 static friction case. Figure 9^ shows the average normal contact 
force. The peak value for the pebble to surface values is due to the base of 
the 'arch' formed by pebbles. The curves for both the 0.65 static friction 
coefficient and the 0.35 static friction coefficient are approximately the same 
because the static friction force is not reaching the full Coulomb limit, so 
both have the same effective /z. 



Figure 9.4 shows the normal times force sums. For the 0.65 case, the peak 
is due to the vibrational impact. For the 0.35 case, the vibration travels deep 
into the reactor bed, producing dust throughout the reactor. Figure [975] shows 
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Figure 9.4: Pebble to Pebble Wear 



the peak dust production coming in the base of the reactor, where the forces 
are the highest, and the greatest lengths are traveled next to the wall. 
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Figure 9.5: Pebble to Surface Wear 



Chapter 10 
Future Work 



The dust production simulation requires both proper dust production wear 
coefficients, and properly determining the correct method of dealing with vi- 
brational issues. It would be useful to determine the number of pebbles that 
need to be simulated to provide a correct representation of a full NGNP-600 
sized reactor. Since the middle portions are geometrically similar, deter- 
mining the amount of recirculation that is required to reach a geometrically 
asymptotic state might allow only a portion of the recirculation to be done. 
Those two changes might allow quicker simulation of full sized reactors. Fi- 
nally, in order to allow sufficiently fast simulations on today's computer hard- 
ware many approximations to the true behavior are done. In the future, some 
of these approximations maybe relaxed. 
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Chapter 11 

Summary and Conclusions 



Research results presented in this dissertation demonstrates a distinct ele- 
ment method that provides high fidelity and yet has reasonable iTin-timcs for 
many pebble fuel element flow simulations. The new static friction test will 
be useful for evaluating any implementation of static friction for spheres. The 
PEBBLES code produced for this dissertation has been able to provide data 
for multiple applications including Dancoff factor calculation, neutronics sim- 
ulation and earthquake simulation. The new static friction model provides 
expected static friction behavior in the reactor including partial matching of 
the Janssen model predictions and correctly matching stability behavior in 
a pyramid. The groundwork has been created for predicting the dust pro- 
duction from wear in a pebble bed reactor once further experimental data 
is available. Future work includes potentially relaxing some of the physical 
approximations made for speed purposes when faster computing hardware 
exists, and investigating new methods for allowing faster simulations. This 
dissertation has provided significant enhancements in simulation of the me- 
chanical movement of pebbles in a pebble bed reactor. 
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Appendix A 

Calculation of Packing 
Fractions 



For determining the volume of a sphere that is inside a vertical slice, the 
following formula can be used 



a = max{—r, hot — z) 
b = min{r, top — z) 

1 



V = Tl 



r\b-a) + -{a^-b^) 
o 



(A.l) 
(A.2) 

(A.3) 



where r is the pebble radius, bot is the bottom of the vertical slice,top is the 
top of the vertical slice and z is the vertical location of the pebble center. 

To determine the area inside a vertical and radial slice, two auxiliary 
functions are defined, one which has the area inside a radial 2d slice, and 
another which has the outside a radial 2d slice. 



Figure A.l shows the area that is inside both a circle of radius c and a 
radial slice of /. The circle is r from the center of the radial circle. Auxiliary 
terms are defined, which include /, the distance from the intersection of the 
segment of the radial circle perpendicularly to the center line, j the distance 
to the intersection of /, 6' the angle of segment, and the angle from the 
segment intersection on the circle side. The area_inside function has the 
following definition: 
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Figure A.l: Area Inside Geometry 




Figure A. 2: Area Outside Geometry 
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Figure A.2 shows the area that is outside the radial shce, but inside the 
circle. The radial slice has a radius of O. The new auxiliary term k is 
the distance from the circle's center to the perpendicular intercept. The 
area_outside function has the following definition: 



Go = 


0.0 iiO>c + r 




(A.13) 


Qo = 


TX(? lie — r > 




(A.14) 


tto = 


T\(? - nO^ if < r H 


- c and < c — r 


(A.15) 


otherwise 
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(A.21) 



Then, the total volume in a radial slice can be determined from the com- 
putation: 



a = 


max{—r, hot — z) 




(A.22) 


b = 


min{r, top — z) 
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Appendix B 

Determination of dust 
production coefficients 



One potential use of the PEBBLES code is to predict the dust production 
of a pebble bed reactor. This section discusses the features that make this 
possible and work that has been done to determine the necessary coefficients. 
Unfortunately, the following literature review discovered a lack of robust wear 
coefficients, which prevents prediction of dust production. 

There are essentially four contact wear mechanisms. Adhesive wear is 
from the contacting surfaces adhesively bonding together, and part of the 
material is pulled away. Abrasive wear is when one of the contacting materials 
is harder than the other, and plows (or shears) away material. Fatigue wear 
is when the surfaces repeatedly contact each other causing fracture of the 
material. The last mechanism is corrosive wear, when chemical corrosion 



causes the surface to behave with increased wear ( Bhushan , 2000 ) . For pebble 
bed reactors, adhesive wear is expected to be the dominate wear mechanism. 
As a first order approximation the adhesive dust production volume is 



(Bhushan 2000): 



N 

V = Kadj^L (B.l) 

In this equation V is the wear volume, Kad is the wear coefficient for 
adhesive wear, L is the length slide and ^ is the real contact area (with 

the normal force and H the hardness). Typically, the hardness and the 
wear coefficient for adhesive wear are combined with the units of either mass 
or volume over force times distance. For two blocks, the length slide is the 
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distance that one of the blocks travels over the other while in contact. Note 
that this formula is only an approximation since the wear volume is only 
approximately linear with respect to both the normal force and the distance 
traveled. Abrasive wear also can be approximated by this model, but fatigue 
and corrosive wear will not be modeled well by this. To the extent that these 
wear mechanisms are present in the pebble bed reactor, this model may be 
less valid. 

The wear coefficient is typically measured by grinding or stroking two 
pieces of graphite against each other, and then measuring the decrease in 
mass. The details of the experiment such as the contact shape and the 
orientation of the relative motion affect the wear coefficient. 

The wear that occurs with graphite depends on multiple factors. A par- 
tial list includes the normal force of contact (load), the temperature of the 
graphite and the past wear history (since wear tends to polish the contact 
surfaces and remove loose grains). The atmosphere that the graphite is in af- 
fects the wear rates since some molecules chemically interact with the carbon 
or are adsorbed on the surface. Neutron damage and other radiation effects 
can damage the structure of the graphite and affect the wear. The type and 
processing of the graphite can affect wear rates. As a related effect, if harder 
and softer graphites interact, the harder one can 'plow' into the softer and 
increase wear rates. 

For graphite on graphite, depending on conditions there can be over three 
orders of magnitude difference in the wear. For example Sheng et al. (2003) 
experimentally determined graphite on graphite in air at room temperature 
can exhibit wear rates of 3.3e-8 g/(Nm) but in the dusting regim^ at 200°C 
the wear coefficient was determined by Lancaster and Pritchard (1980) to 
be 2e-5 g/(Nm) which is about a thousand times greater. For this reason, 
conditions as close to the in-core conditions are needed for determining a 
better approximation of the wear coefficients. 

For tests using nuclear graphite near in-core conditions, the best data 
available to the author is from two independent sets of experiments. One 
data-set emerged from the experiments by Stansfield (1969) and the other 



is from a series of experiments performed at the Tsinghua University (Sheng 



et al. 


2003; 


Luo et al. 


2004, 


2005) 



^In air, above a certain temperature graphite wear transitions to dusting wear, which 
has much greater wear rates. Increased water vapor decreases or ehminates the dusting 
wear. 
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O.M. Stansfield measured friction and wear with different types of graphite 
in hehum at different temperatures (Stansfield, 1969). In the experiments, 
two pieces of graphite were shd against each other hnearly with a 0.32 cm 
stroke. Two different loads were used, one 2-kg mass, and another 8-kg mass. 
The data for wear volumes is only provided graphically, that is, not tabulated, 
therefore only order of magnitude results are available. The wear values were 
about an order of magnitude higher at 25°C than at 400°C and 800°C. There 
was a reduction of friction with increased length slide, but no explanation 
was providec|^ Typical values for the wear rates are lOe-3 cm^/kg for the 
25°C case and lOe-4 cmVkg for the 400°C and 800°C for 12 500 cm distance 
slide. With a density of 1.82 g/cm^, these work out to about 1.5e-6 g/(Nm) 
and 1.5e-7 g/(Nm). These are only about an order of magnitude above room 
temperature wear. 

The second set of experiments were done at the Tsinghua university. 
The first paper measures the wear coefficient of graphite KG-11 via press- 
ing a static specimen against a revolving specimen. The wear is measured 
by weighing the difference in mass before the experiment and after the ex- 
periment. At room temperature in air they measured wear rates of 7.32e-9 
g/ (Nm) with 31 N load with surface contact, 3.29e-8 g/ (Nm) with 31 N load 
with line contact and 3.21e-8 g/(Nm) with 62 N load( |Sheng et al.||2003D . The 
second paper measures the wear coefficient of graphite IG-11 on graphite and 
on steel at varying loads(Luo et al. , 2004). Unfortunately, there are incon- 
sistencies in the units used in the paper. For example, in Table 2 the mean 
wear rate for the lower specimen is listed as 3.0e3 /xg/m, but in the text it is 
listed as 0.3e-3 /ig/m, seven orders of magnitude different. The 30 N of load 
upper specimen wear coefficient for the first 30 minutes is listed as 1.4e-3 
/xg/m, which works out to 4.7e-10 g/(Nm). If 1.4e3 /xg/m is used, this works 
out to 4.7e-4 g/(Nm). Neither of these matches the first paper's results. 
It seems that the units of fig, (or micrograms or l.Oe-6 g) are used where 
mg (or milligrams or l.Oe-3 g) should be. Also, the sign for the exponent 
is inconsistent, where sometimes the negative sign is dropped. These two 
mistakes would make the correct exponents l.Oe-3 mg/m and the measured 
coefficient 1.4e-3 mg/m or 4.7e-7 g/(Nm), which match reasonably well to 
the first paper's values on the order of l.Oe-8 g/(Nm). For the rest of this 
report, it is assumed that these corrections should be used for the Xiaowei 
et al. papers. 



^Possibly this was due to a lubrication effect or the removal of rough or loose surfaces. 
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The third paper measures the temperature effects in hehum(Luo et aL 



2005). The experimental setup is similar to the setup in the second paper, 



but the atmosphere is a helium atmosphere and the temperatures used are 
100°C to 400°C with a load of 30 N. In Figure 2 of that paper, it can be 
qualitatively determined that as the temperature increases, the amount of 
wear increases. As well, the wear tends to have a higher rate initially, and 
then decrease. Since the wear experiment was performed using a 2 mm long 
stroke, it seems plausible that wear rates in an actual pebble bed might be 
closer to the initially higher rates since the pebble flow might be able to 
expose more fresh surfaces of the pebbles to wear. From the graph, there 
does not seem to be a clear trend in the wear as a function of temperature. 
This makes it difficult to estimate wear rates since pebble bed reactor cores 
can have temperatures over 1000° C in normal operation. The highest wear 
rate in Table 2 of the paper is 31.3e-3 mg/m at 30 N, so the highest wear 
rate measured is 1.04e-6 g/(Nm). This is about 20 times lower than wear in 
the dusting regime. Since the total amount of wear (from Fig. 2) between 
200°C and 400°C roughly doubles in the upper specimen and increases by 
approximately 35% in the lower specimen, substantially higher wear rates 
in over 1000°C environments are hard to rule out. Note, however, that the 
opposite temperature trend was observed in the Stansfield paper. 



B.l Calculation of Force in Reactor Bed 

In order to calculate the dust produced in the reactor, the force acting on the 
pebbles is needed. Several different approximations can be used to calculate 
this with varying accuracy. The simplest (but least accurate) method of 
approximating the pressure in the reactor is using the hydrostatic pressure, 
or 



P = pfgh (B.2) 

where P is the pressure at a point, p is the density of the pebbles, / is the 
packing fraction of the pebbles (typical values are near 0.61 or 0.60), g is the 
gravitational acceleration and h is the height below the top of the pebble 
bed. With knowledge of how many contacts there are per unit area or per 
unit volume, this can be converted into pebble to surface or pebble to pebble 
contact forces. This formula is not correct when static friction occurs since 
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the static friction allows forces to be transferred to the walls. Therefore, 



Equation |B.2| over-predicts the actual pressures in the pebble bed. 

In the presence of static friction, more complicated calculations are re- 
quired. The fact that static friction transfers force to the wall was observed 



by the German engineer H.A. Janssen in 1895 (Sperl, 2006). Formulas for the 



pressure on the wall for cylindrical vessels with conical exit chutes were de- 



rived by Walker (1966). Essentially, when the upward force on the wall from 



static friction for a given segment matches the downward gravitational force 
from the additional pebbles in that segment, the pressure stops increasing. 



For a cylinder, the horizontal pressure equation is (Gotoh et al. , 1997): 



1 — exp 



D 



-X 



(B.3) 



where 7 is the bulk weight (or fpg), D is the diameter of the cylinder, //^ 
is the static friction coefficient between the pebbles and the wall, K is the 
Janssen Coefficient, and x is the distance below the top of the pile. 

The Janssen coefficient is dependent upon the pebble to pebble static 
friction coefficient and can be calculated from: 



K 



1 



sm( 



where tan< 



tan ^ n 



sm 



-1 



(B.4) 

1 + sin ^ ' 

Hp and Hp is the pebble to pebble static friction. Since 
, , then K can also be written as: 



K = 2hI-2hpJhI + ^ + ^ 



(B.5) 

The Janssen formula derivations make assumptions that are not necessar- 
ily true for granular materials. These include assuming the granular material 
is a continuum and that the shear forces on the wall are at the Coulomb limit 
(Bratberg et al. , 2005). The static friction force ranges from zero at first con- 



tact up to jj,N (the Coulomb limit) when sufficient shear force has occurred. 
If the force is not at the Coulomb limit, then an effective /x may be able 
to be found and used instead. In general, this assumption will not be true 
when the pebbles are freshly loaded as they will not have slid against the 
wall enough to fully load the static friction. Even after the pebbles have 
been recirculated, they may not reach the Coulomb limit and effective val- 
ues for the static friction constant may be needed instead for predicting the 
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Table B.l: AVR Data 



Name 


Value 


Average Inlet Temperature 


250°C 


Average Outlet iemperature 


950 C 


Pebble Circulation Rate 


300-500 per day 


Dust Produced 


3 kg per year 


Pebbles in Reactor Core 


100,000 


Reactor Radius 


1.5 m 


Outlet Chute Radius 


0.25 m 


Angle of Outlet Cone 


30° 


Control Rod Nose Thickness 


0.3 m 


Radius of Control Rod Nose 


0.15 m 


Feed tube to outlet chute 


2.83 m 



wall pressure. Finally, real reactors have more complicated geometries than 
a smooth cylinder above a cone exit chute. 



B.2 Prior data on dust production 



The 46 MW thermal pebble bed reactor Arbeitsgemeinschaft VersuchsReak- 
tor (AVR) was created in the 1960s in Germany and operated for 21 years. 
The pebbles were added into the reactor through four feeding tubes spaced 
around the reactor and one central feeding tube at the top of the reactor. 
There was one central outlet chute below. Into the reactor cavity there were 
four noses of U shaped graphite with smooth sides for inserting the control 
rods. The cylinder walls contained dimples about 1/2 a pebble diameter 
deep and that alternated location periodically. All the structural graphite 



was a needle coke graphite. Dimensions are shown in Figure B.l and design 
and measured data is provided in Table |B.1[ The measured dust produc- 



tion rate was 3 kg per year. No real conclusions were inferred because of a 
water ingress, an oil ingress, the uncertainty in the composition of the dust 
(i.e., metallic components) and the uncertainty of the location of dust pro- 
duction (Baumer et al. 1990 Atomwirtschaft-Atomtechnik-atw 1966). The 



interior of the AVR reactor reached over 1280° C as determined by melt wire 



experiments( Moormann ) . 



The THTR-300 reactor was a thorium and uranium powered pebble bed 
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Figure B.l: AVR Dimensions 



reactor that first went critical in 1983 and ran through 1988. THTR-300 
produced 16 kg of dust per Full Power Year (FPY), and an estimated 6 kg of 
that was produced in the core of the reactor(Wahsweiler, 1989). The control 
rods in the THTR-300 actually pushed into the pebble bed. On a per pebble 
basis, the amount of dust produced in the THTR-300 is lower than in the 
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Table B.2: THTR Data 



Name 


Value 


Average Inlet Temperature 


250°C 


Average Outlet Temperature 


750°C 


Core Height 


6.0 m 


Pebbles Circulated 


1,300,000 per FPY 


Core Diameter 


5.6 m 


Pebbles in Full Core 


657,000 


Total Dust Produced 


16 kg per FPY 


Estimated In-core Dust 


6 kg per FPY 



AVR. Further data on the THTR-300 is summarized in Table B.2 tht, afb ) 



B.3 Prior Prediction Work 

There are two papers published that attempt to predict the in-core pebble 
dust production. The first paper is "Estimation of Graphite Dust Quantity 



and Size Distribution of Graphite Particle in HTR-10" (Xiaowei et al. 2005) 
and was created to estimate the dust production that the core of the HTR- 
10 reactor would produce. The second is co-authored by this author and 
attempts to estimate the dust that the AVR reactor produced. 

The HTR-10 paper started by calculating from the hydrostatic pressure 
the force between the pebbles at the bottom of the reactor. The force was 
approximated to be SON. The remainder of the paper uses SON as the force 
for conservatism. Note that the HTR-10 paper is in Chinese, so this review 
may contain mistakes in understanding due to language differences. 

The dust production is calculated in three regions, the core of the reactor, 
the outlet chute of the reactor and the fuel loading pipe. As with the other 
papers, the assumption is made that /ig should actually be mg. 

For the core of the reactor the temperature used is 550°C with pebble 
to pebble wear rates of 4.2e-S mg/m extrapolated from 400°C data. The 
pebble to wall wear rates are extrapolated to 480°C to 12.08e-S mg/m from 
the 400° C data. The pebble to pebble wear is estimated to occur foi02.O6 m 
and S.85% of pebbles are estimated to wear against the wall. From this data 



^This is the length shde and is muhipUed by 4.2e-3 mg/m to get per pass dust produc- 
tion 
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the average pebble dust production per pass in the core is determined to be 
8.65e-3 mg for pebble to pebble wear and 0.99e-3 mg from pebble to wall. 
The total in-core graphite dust produced per pebble pass is 9.64e-3 mg. 

The outlet chute wear is estimated to occur for 2.230 m in the graphite 
portion and 1.530 m in the stainless steel portion, and that 44.16% of the 
pebbles wear against the chute. Both these portions are estimated to be at 
400°C. Wear rates of 3.5e-3 mg/m are used for the pebble to pebble wear, 
and 10.4e-3 mg/m for the pebble to graphite chute and 9.7e-3 mg/m for 
pebble to steel. Thus for the outlet chute the upper portion has 18.05e-3 mg 
of dust produced per average pebble and the lower portion has 11.91e-3 mg 
produced for a total outlet chute amount of 29.96e-3 mg. 

The fuel loading pipe is approximately 25 m long and the temperature 
is 200°C which gives a wear value of 2.1e-3 mg/m and 52.50e-3 mg. Thus, 
for an estimated average pebble pass, 10.5% of the dust is produced in core, 
32.5% is produced in the outlet chute and 57.0% is produced in the loading 
pipes. The paper estimates that 50% of the outlet chute graphite dust enters 
the core and that 75% of the graphite dust produced in the fuel loading 
pipes enters the reactor core, for a total amount of graphite dust entering 
the core of 64.0e-3 mg per pebble pass. Since there are 125 pebbles entering 
the reactor a day, and 365 days in a year, this works out to 2.92 g/year of 
pebble dust per year (reported in the paper as 2.74 kg/year due to a precision 
loss and unit errors) ( Xiaowei et al. , 2005 ) . 

HTR-10 has 27 thousand pebbles compared to AVR's 100 thousand and 
a rate of 125 pebbles per day compared to about 400 pebbles per day. A 
crude scaling factor estimate of 35 grams of dust per year would be produced 
per year in AVR. Measured values of dust generation rates from HTR-10 
would provide valuable information on pebble bed reactor dust production 
but appear to be unavailable. 



